Practical Quantitative Finance with R (Automated Stock Trading System)
Table of Contents
- 1. Disclaimer
- 2. Set up
- 3. Examples
- 3.1. Using data table commands
- 3.2. Data Joins
- 3.3. Inner join (would move back to average)
- 3.4. Bollinger Bands
- 3.5. MACD and RSI Indicators
- 3.6. On Balance Volume (OBV) Indicator
- 3.7. Accumulation/Distribution Indicator (A/D)
- 3.8. Williams %R Inidicator
- 3.9. Stochastic Indicator
- 3.10. Commodity Channel Index Indicator (CCI)
- 3.11. Add Custom Indicators
- 4. Chapter 5: Options Pricing
- 5. Chapter 6: Pricing Fixed-Income Instruments
- 6. Chapter 7: Linear Analysis
- 7. Chapter 8: Time Series Analysis
- 8. Chapter 9: Machine Learning
- 9. Trading Strategies and Backtesting
#+Date Created: March 29, 2020 #+Date \today
#+Author website: www.DrXuDotNet.com #+Author contact: jxu@DrXuDotNet.com
1 Disclaimer
1.1 Using Code Examples
You may use the code in this book in your own applications and documentation. You do not need to contact the author for permission unless you are reproducing a significant portion of the code. For example, writing a program that uses several chunks of code from this book does notrequire permission. Selling or distributing the example code listings does require permission. Incorporating a significant amount of example code from this book into your applications and documentation also requires permission. Integrating the example code from this book into commercial products is not allowed without written permission of the author.
2 Set up
2.1 Set working directory and then set it
2.2 Load packages
3 Examples
3.1 Using data table commands
DT<- data.table(iris) head(DT,3) DT[Sepal.Length >=5,mean(Sepal.Width),by=Species] DT[5:8] DT[Species=="setosa"] DT[,mean(Sepal.Length)] #.() command is an alias to list() DT[,.(mean(Sepal.Length),mean(Sepal.Width),mean(Petal.Length),mean(Petal.Width))] DT[,.(avg_Sepal.Length=mean(Sepal.Length),avg_Sepal.Width=mean(Sepal.Width))] DT[Species == "virginica", .(avg_Sepal.Length = mean(Sepal.Length),avg_Sepal.Width = mean(Sepal.Width))] DT[, .(avg_Sepal.Length = mean(Sepal.Length), avg_Sepal.Width = mean(Sepal.Width)), by = Species] DT[Petal.Length>=3, .(avg_Sepal.Length = mean(Sepal.Length), avg_Sepal.Width = mean(Sepal.Width)), by = .(Species, Petal.Width <= 1.5)]
Date FTSE.Open FTSE.High FTSE.Low FTSE.Close FTSE.Rtn N225.Open 1193 2015-01-05 6547.80 6576.74 6404.49 6417.16 -1.9951740 17101.58 1194 2015-01-06 6417.16 6452.66 6328.59 6366.51 -0.7892900 16808.26 1195 2015-01-07 6366.51 6459.74 6366.51 6419.83 0.8375075 17067.40 1196 2015-01-08 6419.83 6580.82 6419.83 6569.96 2.3385354 17318.74 1197 2015-01-09 6569.96 6570.24 6471.38 6501.14 -1.0474950 16970.88 1198 2015-01-13 6501.42 6558.83 6465.19 6542.20 0.6272476 16961.82 1199 2015-01-14 6542.20 6542.20 6353.65 6388.46 -2.3499740 16872.95 1200 2015-01-15 6388.46 6498.78 6298.15 6498.78 1.7268638 16812.96 1201 2015-01-16 6498.78 6553.20 6443.28 6550.27 0.7923026 17000.78 1202 2015-01-19 6550.27 6598.89 6548.00 6585.53 0.5382984 17071.65 1203 2015-01-20 6585.53 6640.44 6585.53 6620.10 0.5249388 17308.72 1204 2015-01-21 6620.10 6728.04 6620.10 6728.04 1.6304890 17306.64 1205 2015-01-22 6728.04 6808.18 6726.24 6796.63 1.0194648 17520.63 1206 2015-01-23 6796.63 6841.73 6796.58 6832.83 0.5326169 17285.71 1207 2015-01-26 6832.83 6855.92 6790.13 6852.40 0.2864113 17649.40 1208 2015-01-27 6852.40 6864.97 6773.54 6811.61 -0.5952659 17615.93 1209 2015-01-28 6811.61 6863.00 6777.08 6825.94 0.2103761 17666.91 1210 2015-01-29 6825.94 6825.94 6750.22 6810.60 -0.2247310 17788.74 1211 2015-01-30 6810.60 6843.98 6749.40 6749.40 -0.8985992 17536.61 1212 2015-02-02 6749.40 6795.52 6731.99 6782.55 0.4911548 17654.60 N225.High N225.Low N225.Close N225.Rtn Direction Predict 1193 17111.36 16881.73 16883.19 -3.01872716 DOWN DOWN 1194 16974.61 16808.26 16885.33 0.01267903 UP DOWN 1195 17243.71 17016.09 17167.10 1.66872385 UP UP 1196 17342.65 17129.53 17197.73 0.17842769 DOWN UP 1197 17087.71 16828.27 17087.71 -0.63973285 UP DOWN 1198 17036.72 16770.56 16795.96 -1.70736737 DOWN DOWN 1199 17141.98 16856.22 17108.70 1.86198505 UP UP 1200 16864.34 16592.57 16864.16 -1.42932586 UP DOWN 1201 17039.80 16911.58 17014.29 0.89022463 UP UP 1202 17366.30 17066.77 17366.30 2.06891817 UP UP 1203 17329.03 17181.55 17280.48 -0.49417728 UP DOWN 1204 17355.74 17229.21 17329.02 0.28088954 UP UP 1205 17532.06 17460.76 17511.75 1.05447668 UP UP 1206 17471.94 17285.71 17468.52 -0.24686550 UP DOWN 1207 17768.41 17633.47 17768.30 1.71612282 DOWN UP 1208 17850.59 17615.93 17795.73 0.15437429 UP DOWN 1209 17778.83 17575.10 17606.22 -1.06491704 DOWN DOWN 1210 17808.47 17661.10 17674.39 0.38719225 DOWN DOWN 1211 17628.40 17459.45 17558.04 -0.65830593 UP DOWN 1212 17654.60 17271.87 17335.85 -1.26545711 UP DOWN
3.2 Data Joins
symbols<-data.table(id=c(1:4),ticker=c("AMAT","IBM","INTC","MSFT")) prices<-data.table(id=c(1,1,1,2,2,2,3,3,3,5,5,5), date = rep(c('2016-01-04','2016-01-05','2016-01-06'),4), price=c(18.47,18.49,17.73,135.95,135.85,135.17,33.99,33.83,33.08,54.80,55.05,54.05)) setkey(symbols,id) setkey(prices,id)
1 | 2016-01-04 | 18.47 |
1 | 2016-01-05 | 18.49 |
1 | 2016-01-06 | 17.73 |
2 | 2016-01-04 | 135.95 |
2 | 2016-01-05 | 135.85 |
2 | 2016-01-06 | 135.17 |
3 | 2016-01-04 | 33.99 |
3 | 2016-01-05 | 33.83 |
3 | 2016-01-06 | 33.08 |
5 | 2016-01-04 | 54.8 |
5 | 2016-01-05 | 55.05 |
5 | 2016-01-06 | 54.05 |
3.3 Inner join (would move back to average)
3.4 Bollinger Bands
gs<-getSymbols("GS",from="2019-06-01",auto.assign=FALSE) chartSeries(gs,theme="white",minor.ticks=FALSE,TA = "addEnvelope(20,5,'SMA')") chartSeries(gs,theme="white",minor.ticks=FALSE,TA = "addBBands(n=20,sd=2,maType='EMA')")
3.5 MACD and RSI Indicators
3.5.1 MACD Description
Measures difference between short and long EMAs plotted alongside short 9-day EMA of the MACD Trading Strategy: Look at which side of zero tthe MACD lines are on. Above 0 for a long time–trend is likely going up; below 0, trend likely going down Potential buy when MACD moves above 0; Potential sell when MACD moves below 0
3.5.2 Relative Strength Index Description
RSI is a momentum oscillator, calculating velocity and absolute stock price strength. Failure swing occurs when RSI enters above (or below) 70 (or 30) and then crosses back through that line.
chartSeries(gs,theme="white",minor.ticks=FALSE,TA = "addMACD(); addRSI(); addVo()")
3.6 On Balance Volume (OBV) Indicator
3.6.1 Description
OBV shows momentum of volume that can be used to predict stock price changes.
"When volume increases sharply without a significant change in teh stock price, the price will eventually jump upward, and vice versa."
3.6.2 Theory
Basically OBV change come before price changes because when public sees 'smart money' flowing in, they also buy shares to flow in. If stock prices changes precede OBV changes, non-confirmations occurred, which can occur when stock rises without or before the OBV (bulle market top) or when the stock falls without or before the OBV (bear market bottoms. If OBV changes to increasing or decreasing trends, OBV breakouts occur, and OBV breakouts usually precede price breakouts. But long on upside OBV breakouts and sell short on downside OBV breakouts.
chartSeries(gs,theme="white",minor.ticks=FALSE,TA = "addOBV(); addVo()")
3.7 Accumulation/Distribution Indicator (A/D)
3.7.1 Description
Volume changes weight price changes. This is a variant of OBV. If A/D increases, then many people are accumulating stock; decreaes mean distibuting stock. Different values for A/D and price indicate upcoming price changes. Usually, price changes move in same direction and A/D changes. That is, price follows direction of A/D.
#+name AD
chartSeries(gs,theme="white",minor.ticks=FALSE,TA = "addChAD(); addVo()")
3.8 Williams %R Inidicator
3.8.1 Description
Momentum indicator used to measure overbought or oversold. Sell when %R reaches -10% or highest and buy when it reaches -90% or lower. Works best in trending markets (bear or bull).
chartSeries(gs,theme="white",minor.ticks=FALSE,TA = "addWPR(n=10); addVo()")
3.9 Stochastic Indicator
3.9.1 Description
Smooths a simple stochastic indicator as a moving everage of a simple stochastic indicator. The indicator tries to preict price turning points by comparing the stock close price to its price range.
3.9.2 Theory
In upwward-trending markets, prices close near their high and in downward-trending markets prices close near their low. Thus, transaction signals occur when %K (simple stochastic) cross through %D (moving average).
3.9.3 Three popular trading rules
- Buy when indicator falls below specific level and then rises above (e.g., 20%). Sell indicator rises above specific level then falls below (80%).
- Buy when %K line rises above the %D line. Sell when %K line fallse below the %D line.
- Divergences occur when priceses are making a series of new high and the indicator is failing to surpass its previous highs and vice versa.
chartSeries(gs,name="GS",theme="white",minor.ticks=FALSE,TA = "addVo()") addTA(SMI(HLC(gs)),col=c("blue","red"),lwd=c(2,2))
3.10 Commodity Channel Index Indicator (CCI)
3.10.1 Description:
Mean deviation of daily average price from moving average. A value > 100 indicates overbought and < -100 indicates oversold. Used by traders to identify cyclical trends.
3.10.2 Theory
Stock prices move in cycles. Buy when > 100 and closed if move back below 100. Sell when < -100 and then close this position when CCI moves back above -100.
#+name cci
chartSeries(gs,name="GS",theme="white",minor.ticks=FALSE,TA = "addVo()") chartSeries(gs,name="GS",theme="white",minor.ticks=FALSE,TA = "addCCI(15); addVo()")
3.11 Add Custom Indicators
Dr. Xu provides code for a moving linear regression close price, added to the plotHelper.R script.
rr<-selfLinearRegression(gs,20) slr<-rr[,1] + rr[,2] * as.numeric(index(gs)) slrDiff<-gs[,4] - slr chartSeries(gs,name="GS",theme="white",TA=NULL) addTA(slr,on=1,col="blue",lwd=2) addTA(slrDiff,col="red",lwd=2)
After adding them to plothelper file, i will now convert these to generic TA functions
addSLR<-newTA(FUN=SLR,col="blue",lwd=2,on=1) addSLRDiff<-newTA(FUN=SLRDiff,col="red",lwd=2) chartSeries(gs,name="GS",theme="white",minor.ticks=FALSE,TA="addSLR(5);addSLRDiff(5)") amrn <- getSymbols("AMRN",from="2019-01-01",auto.assign=FALSE) chartSeries(amrn,name="AMRN",theme="white",minor.ticks=FALSE,TA="addSLR(5);addSLRDiff(5)")
4 Chapter 5: Options Pricing
4.1 Pricing European Options
4.1.1 Black-Scholes Model
Run black scholes model function using default parameter values
call = blackScholes(optionType="c") call put<-blackScholes(optionType="p") put
[1] 12.95233 [1] 9.259623
[1] 12.95233 [1] 9.259623
# can also use vector maturities<-seq(.1,1,.1) call<-blackScholes(optionType="c",maturity=maturities) put<-blackScholes(optionType="p",maturity=maturities) df<-data.frame(Maturity=maturities,Call=call,Put=put);df
4.1.2 Black-Scholes Greeks
- Delta
Measures rate of change of option price relative to underlying asset's price (partial derivative). Also provides the probability an option will expire in-the-money (ITM). One way to measure potential returns on a trade is to compare delta to option prices across different strike prices.
maturities<-seq(.1,1,.1) callDelta<-blackScholesDelta(optionType="c",maturity=maturities) putDelta<-blackScholesDelta(optionType = "p",maturity=maturities) delta<-data.frame(Maturity = maturities,CallDelta=callDelta,PutDelta=putDelta) delta
- Gamma:
Rate of change of Delta to small changes in the underlying asset's price
maturities<-seq(.1,1,.1) gamma<-blackScholesGamma(maturity=maturities) # uses default input paramters gamma<-data.frame(Maturity=maturities,Gamma=gamma) gamma
- Theta
Rate of change of option price relative to small change in time (time decays)
maturities<-seq(.1,1,.1) ctheta<-blackScholesTheta(optionType="c",maturity=maturities) ptheta<-blackScholesTheta(optionType="p",maturity=maturities) theta<-data.frame(Maturity=maturities,CallTheta=ctheta,PutTheta=ptheta) theta
- Rho
Rate of change of option price rlative to unit change in risk-free interest rate
- Description
Rho is a risk measure that simply tells us by how much call and put prices change as a result of changing interest rates. Rho for ITM option will be largest because of arbitrary activity. (?)
maturities<-seq(.1,1,.1) crho<-blackScholesRho(optionType="c",maturity=maturities) prho<-blackScholesRho(optionType="p",maturity=maturities) rho<-data.frame(Maturity=maturities,CallRho=crho,PutRho=prho) rho
- Description
- Vega
Rate of change of option price relative to chagne in volatitlity of underlying asset Measure of Risk
maturities<-seq(.1,1,.1) vega<-blackScholesVega(maturity=maturities) vega<-data.frame(Maturity=maturities,Vega = vega) vega
- Implied Volatility [FUNCTION NEEDS TO BE VECTORIZED]
Uses bisection method to calculate implied volatility, which is used when Vega us unknown analytically (as is the case with American Options) I think this function is not vectorized??
vol<-volTest() # He used default values for this test, but it seems like he created the skeleton to vectorize this vol
- 3D Plot of implied volatility [NEEDS TO BE VECTORIZED]
# define arbitrary vectors x<-seq(.1,1,length=30) y<-seq(9.5,10.5,length=30) # create 3D graph for call option cvol <-volTest3d(optionType = "c", x, y) persp3Drgl(x = x, y = y, z = cvol, border="black", bty = "b2", ticktype="detailed", xlab = "Maturity", ylab = "Strike", zlab = "cVol") plotdev() # this function appears necessary!!
# create 3D graph for put option pvol<-volTest3d(optionType="p",x,y) persp3Drgl(x=x,y=y,z=pvol,border="black",bty="b2",ticktype="detailed", xlab="Maturity",ylab="Strike",zlab="pVol") plotdev()
- 3D Plot of implied volatility [NEEDS TO BE VECTORIZED]
4.2 Pricing American Options
4.2.1 Note
No analytical solutions, only closed-form approximations. This book uses Barone-Adesi and Whaley (BAW) method.
BAW approximation. [may need to be vectorized] This serves as a direct comparison of this calculation and with results from the literature BAW value is from the literature (which study?); calculation is our function
df<-americanOptionTest() df
4.2.2 Pricing Barrier Options
Here, we will discuss how to price the standard barrier options using the formulas provided in the literature (D. R. Rich, Advances in Futures and Options Research, 7, 267, 1994, and E. G. Haug's book The complete Guide to Option Pricing Formulas).
# [may need to be vectorized] barrierOptionTest("c") barrierOptionTest("p") barrierOptionTest("c",barrierLevel=110) barrierOptionTest("p",barrierLevel=110)
4.2.3 Pricing European Options Using RQuantLib
The main reason to use QuantLib in this book is that QuantLib provides not only the option-pricing tool, but also comprehensive frameworks for pricing fixed-income instruments, quantitative analysis, modeling, trading, and risk managementof financial assets. The books suggests creating a SWIG (simplified Wrapper and Interface Generator) R version of QuantLib, because the actual QuantLib has many many more features.
RQuantLib is not compatibile with R v3.6.0, however, I think this code is important so I will copy it in
- Scalar inputs
library(RQuantLib) eo <-EuropeanOption(type="call", underlying = 100, strike = 100, dividendYield =0.06, riskFreeRate = 0.1, maturity = 0.5, volatility = 0.3) eo
OptionType Strike Spot DivYield Rate Maturity Vol BAWValue Calculated 1 Call 100 90 0.1 0.1 0.1 0.15 0.0206 0.02063525 2 Call 100 100 0.1 0.1 0.1 0.15 1.8771 1.87693090 3 Call 100 110 0.1 0.1 0.1 0.15 10.0089 10.00607824 4 Call 100 90 0.1 0.1 0.1 0.25 0.3159 0.31590056 5 Call 100 100 0.1 0.1 0.1 0.25 3.1280 3.12769717 6 Call 100 110 0.1 0.1 0.1 0.25 10.3919 10.39013102 7 Call 100 90 0.1 0.1 0.1 0.35 0.9495 0.94946687 8 Call 100 100 0.1 0.1 0.1 0.35 4.3777 4.37768505 9 Call 100 110 0.1 0.1 0.1 0.35 11.1679 11.16785412 10 Call 100 90 0.1 0.1 0.5 0.15 0.8208 0.82081816 11 Call 100 100 0.1 0.1 0.5 0.15 4.0842 4.08413110 12 Call 100 110 0.1 0.1 0.5 0.15 10.8087 10.80852656 13 Call 100 90 0.1 0.1 0.5 0.25 2.7437 2.74360700 14 Call 100 100 0.1 0.1 0.5 0.25 6.8015 6.80135535 15 Call 100 110 0.1 0.1 0.5 0.25 13.0170 13.01674141 16 Call 100 90 0.1 0.1 0.5 0.35 5.0063 5.00616776 17 Call 100 100 0.1 0.1 0.5 0.35 9.5106 9.51031688 18 Call 100 110 0.1 0.1 0.5 0.35 15.5689 15.56843626 19 Put 100 90 0.1 0.1 0.1 0.15 10.0000 10.00000000 20 Put 100 100 0.1 0.1 0.1 0.15 1.8770 1.87693172 21 Put 100 110 0.1 0.1 0.1 0.15 0.0410 0.04099685 22 Put 100 90 0.1 0.1 0.1 0.25 10.2533 10.25300802 23 Put 100 100 0.1 0.1 0.1 0.25 3.1277 3.12769785 24 Put 100 110 0.1 0.1 0.1 0.25 0.4562 0.45618386 25 Put 100 90 0.1 0.1 0.1 0.35 10.8787 10.87850494 26 Put 100 100 0.1 0.1 0.1 0.35 4.3777 4.37768536 27 Put 100 110 0.1 0.1 0.1 0.35 1.2402 1.24020578 28 Put 100 90 0.1 0.1 0.5 0.15 10.5595 10.55920803 29 Put 100 100 0.1 0.1 0.5 0.15 4.0842 4.08412981 30 Put 100 110 0.1 0.1 0.5 0.15 1.0822 1.08219870 31 Put 100 90 0.1 0.1 0.5 0.25 12.4419 12.44164705 32 Put 100 100 0.1 0.1 0.5 0.25 6.8014 6.80135335 33 Put 100 110 0.1 0.1 0.5 0.25 3.3226 3.32258007 34 Put 100 90 0.1 0.1 0.5 0.35 14.6945 14.69432660 35 Put 100 100 0.1 0.1 0.5 0.35 9.5104 9.51031428 36 Put 100 110 0.1 0.1 0.5 0.35 5.8823 5.88219585 Maturity DownIn DownOut UpIn UpOut 1 0.1 0.04440666 3.910483 NA NA 2 0.2 0.34306327 5.323445 NA NA 3 0.3 0.80206291 6.193476 NA NA 4 0.4 1.32196272 6.799050 NA NA 5 0.5 1.86068640 7.252225 NA NA 6 0.6 2.39982948 7.607623 NA NA 7 0.7 2.93086961 7.895483 NA NA 8 0.8 3.44977491 8.134120 NA NA 9 0.9 3.95470307 8.335396 NA NA 10 1.0 4.44493474 8.507400 NA NA Maturity DownIn DownOut UpIn UpOut 1 0.1 2.658598 0.89947884 NA NA 2 0.2 4.435216 0.44398863 NA NA 3 0.3 5.553516 0.27047255 NA NA 4 0.4 6.386166 0.18522032 NA NA 5 0.5 7.054973 0.13632790 NA NA 6 0.6 7.614533 0.10534363 NA NA 7 0.7 8.094459 0.08429773 NA NA 8 0.8 8.512895 0.06925608 NA NA 9 0.9 8.881928 0.05807957 NA NA 10 1.0 9.210106 0.04951676 NA NA Maturity DownIn DownOut UpIn UpOut 1 0.1 NA NA 3.243091 0.71179838 2 0.2 NA NA 5.334045 0.33246351 3 0.3 NA NA 6.797205 0.19833355 4 0.4 NA NA 7.986674 0.13433890 5 0.5 NA NA 9.014690 0.09822134 6 0.6 NA NA 9.931888 0.07556512 7 0.7 NA NA 10.766068 0.06028414 8 0.8 NA NA 11.534477 0.04941784 9 0.9 NA NA 12.248726 0.04137290 10 1.0 NA NA 12.917110 0.03522427 Maturity DownIn DownOut UpIn UpOut 1 0.1 NA NA 0.07561793 3.482459 2 0.2 NA NA 0.43960224 4.439602 3 0.3 NA NA 0.91278267 4.911206 4 0.4 NA NA 1.39847905 5.172907 5 0.5 NA NA 1.86676853 5.324532 6 0.6 NA NA 2.30858479 5.411292 7 0.7 NA NA 2.72198629 5.456770 8 0.8 NA NA 3.10762902 5.474522 9 0.9 NA NA 3.46709922 5.472908 10 1.0 NA NA 3.80225337 5.457369 Concise summary of valuation for EuropeanOption value delta gamma vega theta rho divRho 9.1129 0.5623 0.0179 26.8318 -9.3873 23.5571 -28.1136
- Vector Inputs
spot = seq(10,200,length=30) mat = seq(0.1,3,length=30) call <-EuropeanOptionArrays(type = "call",strike = 100, dividendYield = 0.06, riskFreeRate = 0.1, volatility = 0.3, underlying = spot, maturity = mat) put <-EuropeanOptionArrays(type = "put",strike = 100, dividendYield = 0.06, riskFreeRate = 0.1, volatility = 0.3, underlying = spot, maturity = mat) call$parameters.type call$parameters.strike call$parameters.underlying
[1] "call" [1] 100 [1] 10.00000 16.55172 23.10345 29.65517 36.20690 42.75862 49.31034 [8] 55.86207 62.41379 68.96552 75.51724 82.06897 88.62069 95.17241 [15] 101.72414 108.27586 114.82759 121.37931 127.93103 134.48276 141.03448 [22] 147.58621 154.13793 160.68966 167.24138 173.79310 180.34483 186.89655 [29] 193.44828 200.00000
- Option Prices
head(call$value,2) persp3Drgl(x= spot, y = mat, z = call$value, border="black", bty = "b2", ticktype = "detailed", xlab = "spot", ylab ="maturity", zlab = "cPrice") plotdev()
persp3Drgl(x= spot, y = mat, z = put$value, border="black", bty = "b2", ticktype = "detailed", xlab = "spot", ylab = "maturity", zlab = "pPrice") plotdev()
- Delta
persp3Drgl(x= spot, y = mat, z = call$delta, border="black", bty = "b2", ticktype = "detailed", xlab = "spot", ylab = "maturity", zlab = "cDelta") plotdev()
persp3Drgl(x=spot, y = mat, z = put$delta, border="black", bty = "b2", ticktype = "detailed", xlab = "spot", ylab = "maturity", zlab = "pDelta") plotdev()
image2D(z=call$delta,x=spot,y=mat, contour = T, alpha=0.7, resfac = 5, xlab ="spot", ylab = "maturity")
persp3Drgl(x= spot, y = mat, z = call$gamma, border="black", bty = "b2", ticktype= "detailed", xlab = "spot", ylab = "maturity", zlab = "Gamma") plotdev()
image2D(z=call$gamma,x=spot,y=mat, contour = list(nlevels = 20), alpha=0.7,resfac = 5, xlab ="spot", ylab = "maturity")
- Theta
persp3Drgl(x= spot, y = mat, z = call$theta, border="black", bty = "b2", ticktype= "detailed", xlab = "spot", ylab = "maturity", zlab = "cTheta") plotdev()
persp3Drgl(x= spot, y = mat, z = put$theta, border="black", bty = "b2", ticktype = "detailed", xlab = "spot", ylab = "maturity", zlab = "pTheta") plotdev()
image2D(z=call$theta,x = spot,y = mat, contour = list(nlevels = 20), alpha=0.7, resfac = 5, xlab = "spot", ylab = "maturity")
image2D(z = put$theta,x = spot,y = mat, contour = list(nlevels = 20), alpha=0.7,resfac = 5, xlab = "spot", ylab = "maturity")
- Vega
persp3Drgl(x= spot, y = mat, z = put$vega, border="black", bty = "b2", ticktype= "detailed", xlab = "spot", ylab = "maturity", zlab = "Vega") plotdev()
image2D(z = put$vega,x = spot,y =mat, contour = list(nlevels = 10), alpha=0.7, resfac = 5, xlab = "spot", ylab = "maturity")
persp3Drgl(x= spot, y = mat, z = call$rho, border="black", bty = "b2", ticktype = "detailed", xlab = "spot", ylab = "maturity", zlab = "cRho") plotdev()
persp3Drgl(x= spot, y = mat, z = put$rho, border="black", bty = "b2", ticktype= "detailed", xlab = "spot", ylab = "maturity", zlab = "pRho") plotdev()
image2D(z = call$rho,x = spot,y = mat, contour = list(nlevels = 10), alpha=0.7, resfac = 5, xlab = "spot", ylab = "maturity")
image2D(z = put$rho,x = spot,y = mat, contour = list(nlevels = 10), alpha=0.7, resfac = 5, xlab = "spot", ylab = "maturity")
persp3Drgl(x = spot, y = mat, z = call$divRho, border="black", bty = "b2", ticktype = "detailed", xlab = "spot", ylab = "maturity", zlab = "cDivRho") plotdev()
persp3Drgl(x = spot, y = mat, z = put$divRho, border="black", bty = "b2", ticktype = "detailed", xlab = "spot", ylab = "maturity", zlab = "pDivRho") plotdev()
image2D(z = call$divRho,x = spot,y = mat, contour = list(nlevels = 10), alpha=0.7, resfac = 5, xlab = "spot", ylab = "maturity")
image2D(z = put$divRho,x = spot,y = mat, contour = list(nlevels = 10), alpha=0.7, resfac = 5, xlab = "spot", ylab = "maturity")
- Implied Volatility
volTestRQuantLib()
- Pricing American Options Using RQuantLib
americanOptionTestRQuantLib() # americanoption function is part of RQuantLib package
- Greeks
americanOptionGreeks()
- Implied Volatility
americanOptionvolTest()
- Pricing Barrier Options Using RQuantLib
barrierOptionTestRQuantLib(optionType="c",barrierLevel=90) barrierOptionTestRQuantLib(optionType="p",barrierLevel=90) barrierOptionTestRQuantLib(optionType="c",barrierLevel=110) barrierOptionTestRQuantLib(optionType="p",barrierLevel=110)
If you want to value options in real-world market with discrete dividend payments using the interest rate curve and volatility surface, you will need to bootstrap the required dividend, interest rate, and volatility curves from real market data.
Since the current version of RQuantLib does not exposes the interface to access necessary features, you need to use the original QuantLib or its equivalent SWIG R version to evaluate real-world options.
If you are interested in this topic, my previous published book "Practical C# and WPF for Financial Markets" provided an example on how to price a real-world American option using the C# SWIG version of QuantLib.
This example will give you an idea for what is required to price options in the real-world market.
By following the steps and templates described in that example, you should be able to price various options in real-world markets.
5 Chapter 6: Pricing Fixed-Income Instruments
5.1 Simple Bonds Pricing
5.2 Pricing Bonds Using RQuantLib (p. 185, pdf 205)
source("./R/fixedIncomeHelper.R")
# parameters: evalDate, issuedate, maturity, settlementDays, faceValue, rate, coupon, frequency) library(RQuantLib) library(dplyr) bondFlatRate("2015-12-16", "2015-12-16", "2025-12-16", 1,1000, 0.06,0.05,"Annual")
Maturity Price CallVol PutVol 1 0.1 0.15 0.2508532 NA 2 0.2 0.20 0.1971808 NA 3 0.3 0.25 0.1756172 NA 4 0.4 0.30 0.1638570 0.05729459 5 0.5 0.35 0.1565258 0.08812524 6 0.6 0.40 0.1516229 0.10983901 7 0.7 0.45 0.1482221 0.12780934 8 0.8 0.50 0.1458320 0.14361290 9 0.9 0.55 0.1441666 0.15797678 10 1.0 0.60 0.1430476 0.17130630 Type Strike Spot Yield Rate Maturity Vol BAW BAWCalc CN 1 call 100 90 0.1 0.1 0.1 0.15 0.0206 0.02063560 0.02051471 2 call 100 100 0.1 0.1 0.1 0.15 1.8771 1.87692069 1.87637235 3 call 100 110 0.1 0.1 0.1 0.15 10.0089 10.00606009 10.01038104 4 call 100 90 0.1 0.1 0.1 0.25 0.3159 0.31590017 0.31522682 5 call 100 100 0.1 0.1 0.1 0.25 3.1280 3.12768316 3.12676883 6 call 100 110 0.1 0.1 0.1 0.25 10.3919 10.39012353 10.39647379 7 call 100 90 0.1 0.1 0.1 0.35 0.9495 0.94947033 0.94813190 8 call 100 100 0.1 0.1 0.1 0.35 4.3777 4.37766904 4.37638710 9 call 100 110 0.1 0.1 0.1 0.35 11.1679 11.16785000 11.17235572 10 call 100 90 0.1 0.1 0.5 0.15 0.8208 0.82082169 0.81125463 11 call 100 100 0.1 0.1 0.5 0.15 4.0842 4.08411689 4.06806946 12 call 100 110 0.1 0.1 0.5 0.15 10.8087 10.80852655 10.80993671 13 call 100 90 0.1 0.1 0.5 0.25 2.7437 2.74360046 2.72202909 14 call 100 100 0.1 0.1 0.5 0.25 6.8015 6.80134160 6.77456288 15 call 100 110 0.1 0.1 0.5 0.25 13.0170 13.01673013 13.00035432 16 call 100 90 0.1 0.1 0.5 0.35 5.0063 5.00615862 4.97310988 17 call 100 100 0.1 0.1 0.5 0.35 9.5106 9.51030764 9.47274597 18 call 100 110 0.1 0.1 0.5 0.35 15.5689 15.56842765 15.53773254 19 put 100 90 0.1 0.1 0.1 0.15 10.0000 10.00000000 10.00037680 20 put 100 100 0.1 0.1 0.1 0.15 1.8770 1.87692247 1.87636987 21 put 100 110 0.1 0.1 0.1 0.15 0.0410 0.04099631 0.04081541 22 put 100 90 0.1 0.1 0.1 0.25 10.2533 10.25300774 10.25977819 23 put 100 100 0.1 0.1 0.1 0.25 3.1277 3.12768538 3.12676196 24 put 100 110 0.1 0.1 0.1 0.25 0.4562 0.45618539 0.45535459 25 put 100 90 0.1 0.1 0.1 0.35 10.8787 10.87850882 10.88351168 26 put 100 100 0.1 0.1 0.1 0.35 4.3777 4.37767141 4.37637377 27 put 100 110 0.1 0.1 0.1 0.35 1.2402 1.24020907 1.23866797 28 put 100 90 0.1 0.1 0.5 0.15 10.5595 10.55921289 10.56350493 29 put 100 100 0.1 0.1 0.5 0.15 4.0842 4.08411682 4.06806423 30 put 100 110 0.1 0.1 0.5 0.15 1.0822 1.08220238 1.07075143 31 put 100 90 0.1 0.1 0.5 0.25 12.4419 12.44164207 12.42895237 32 put 100 100 0.1 0.1 0.5 0.25 6.8014 6.80134134 6.77454946 33 put 100 110 0.1 0.1 0.5 0.25 3.3226 3.32257190 3.29798979 34 put 100 90 0.1 0.1 0.5 0.35 14.6945 14.69431918 14.66800414 35 put 100 100 0.1 0.1 0.5 0.35 9.5104 9.51030684 9.47272287 36 put 100 110 0.1 0.1 0.5 0.35 5.8823 5.88218976 5.84511589 Type Strike Spot Yield Rate Maturity Vol Delta Gamma 1 call 100 90 0.1 0.1 0.1 0.15 0.01385890 0.008277720 2 call 100 100 0.1 0.1 0.1 0.15 0.50568750 0.083708772 3 call 100 110 0.1 0.1 0.1 0.15 0.98493010 0.012292632 4 call 100 90 0.1 0.1 0.1 0.25 0.09706771 0.024087949 5 call 100 100 0.1 0.1 0.1 0.25 0.51194080 0.050200312 6 call 100 110 0.1 0.1 0.1 0.25 0.89147570 0.021817980 7 call 100 90 0.1 0.1 0.1 0.35 0.18332460 0.026575630 8 call 100 100 0.1 0.1 0.1 0.35 0.51819042 0.035830658 9 call 100 110 0.1 0.1 0.1 0.35 0.81669675 0.021785914 10 call 100 90 0.1 0.1 0.5 0.15 0.16627961 0.025843144 11 call 100 100 0.1 0.1 0.5 0.15 0.50442954 0.037112940 12 call 100 110 0.1 0.1 0.5 0.15 0.82146404 0.024228163 13 call 100 90 0.1 0.1 0.5 0.25 0.29392771 0.021344916 14 call 100 100 0.1 0.1 0.5 0.25 0.51797280 0.022214635 15 call 100 110 0.1 0.1 0.5 0.25 0.71913541 0.017470905 16 call 100 90 0.1 0.1 0.5 0.35 0.36702944 0.016638068 17 call 100 100 0.1 0.1 0.5 0.35 0.53147895 0.015810855 18 call 100 110 0.1 0.1 0.5 0.35 0.67696829 0.013105187 19 put 100 90 0.1 0.1 0.1 0.15 -0.99624693 0.017040768 20 put 100 100 0.1 0.1 0.1 0.15 -0.48692115 0.083708893 21 put 100 110 0.1 0.1 0.1 0.15 -0.02082302 0.009589597 22 put 100 90 0.1 0.1 0.1 0.25 -0.90139918 0.025405092 23 put 100 100 0.1 0.1 0.1 0.25 -0.48066874 0.050200435 24 put 100 110 0.1 0.1 0.1 0.25 -0.10557336 0.020950234 25 put 100 90 0.1 0.1 0.1 0.35 -0.81220420 0.027218597 26 put 100 100 0.1 0.1 0.1 0.35 -0.47442045 0.035830782 27 put 100 110 0.1 0.1 0.1 0.35 -0.17813424 0.021359560 28 put 100 90 0.1 0.1 0.5 0.15 -0.82418316 0.030816396 29 put 100 100 0.1 0.1 0.5 0.15 -0.46374332 0.037113046 30 put 100 110 0.1 0.1 0.5 0.15 -0.16363443 0.020944568 31 put 100 90 0.1 0.1 0.5 0.25 -0.68209262 0.023051171 32 put 100 100 0.1 0.1 0.5 0.25 -0.45021800 0.022214748 33 put 100 110 0.1 0.1 0.5 0.25 -0.25502885 0.016335028 34 put 100 90 0.1 0.1 0.5 0.35 -0.60516280 0.017498071 35 put 100 100 0.1 0.1 0.5 0.35 -0.43673870 0.015810977 36 put 100 110 0.1 0.1 0.5 0.35 -0.29428567 0.012531174 Maturity Price CallVol PutVol 1 0.1 0.15 0.2508051 NA 2 0.2 0.20 0.1971306 NA 3 0.3 0.25 0.1755628 NA 4 0.4 0.30 0.1637980 NA 5 0.5 0.35 0.1564623 NA 6 0.6 0.40 0.1515546 NA 7 0.7 0.45 0.1481491 NA 8 0.8 0.50 0.1457541 0.06362939 9 0.9 0.55 0.1440837 0.11678677 10 1.0 0.60 0.1429595 0.13284692 11 1.1 0.65 0.1422649 0.14625304 12 1.2 0.70 0.1419202 0.15825871 13 1.3 0.75 0.1418685 0.16929832 14 1.4 0.80 0.1420697 0.17962205 15 1.5 0.85 0.1424929 0.18939537 Maturity DownIn DownOut UpIn UpOut 1 0.1 0.0444071 3.910475 NA NA 2 0.2 0.3430603 5.323442 NA NA 3 0.3 0.8020655 6.193470 NA NA 4 0.4 1.3219689 6.799042 NA NA 5 0.5 1.8606924 7.252219 NA NA 6 0.6 2.3998328 7.607621 NA NA 7 0.7 2.9308691 7.895485 NA NA 8 0.8 3.4497707 8.134126 NA NA 9 0.9 3.9546957 8.335406 NA NA 10 1.0 4.4449251 8.507412 NA NA Maturity DownIn DownOut UpIn UpOut 1 0.1 2.658597 0.89947200 NA NA 2 0.2 4.435218 0.44398095 NA NA 3 0.3 5.553525 0.27045990 NA NA 4 0.4 6.386176 0.18520844 NA NA 5 0.5 7.054978 0.13632256 NA NA 6 0.6 7.614530 0.10534789 NA NA 7 0.7 8.094446 0.08431218 NA NA 8 0.8 8.512873 0.06927962 NA NA 9 0.9 8.881899 0.05811001 NA NA 10 1.0 9.210074 0.04955122 NA NA Maturity DownIn DownOut UpIn UpOut 1 0.1 NA NA 3.243090 0.71179159 2 0.2 NA NA 5.334038 0.33246421 3 0.3 NA NA 6.797204 0.19833079 4 0.4 NA NA 7.986681 0.13432968 5 0.5 NA NA 9.014705 0.09820570 6 0.6 NA NA 9.931909 0.07554411 7 0.7 NA NA 10.766095 0.06025935 8 0.8 NA NA 11.534506 0.04939127 9 0.9 NA NA 12.248755 0.04134689 10 1.0 NA NA 12.917135 0.03520146 Maturity DownIn DownOut UpIn UpOut 1 0.1 NA NA 0.07561724 3.482451 2 0.2 NA NA 0.43959907 4.439600 3 0.3 NA NA 0.91278507 4.911200 4 0.4 NA NA 1.39848564 5.172898 5 0.5 NA NA 1.86677683 5.324523 6 0.6 NA NA 2.30859306 5.411284 7 0.7 NA NA 2.72199355 5.456764 8 0.8 NA NA 3.10763477 5.474518 9 0.9 NA NA 3.46710326 5.472906 10 1.0 NA NA 3.80225568 5.457370 Error in Date(as.numeric(format(rDate, "%d")), format(rDate, "%B"), as.numeric(format(rDate, (from fixedIncomeHelper.R#168) : could not find function "Date"
5.3 Pricing Bonds Using QuantLib-SWIG R
bondFlatRate("2015-12-16", "2015-12-16", "2025-12-16", 1, 1000,0.05,0.05,"Annual")
Error in Date(as.numeric(format(rDate, "%d")), format(rDate, "%B"), as.numeric(format(rDate, : could not find function "Date"
5.3.1 Pricing Bonds with a Rate Curve
#+call fixedIncomeHelpery
bondCurveRate()
Error in Date(15, "January", 2015) : could not find function "Date"
5.3.2 Zero-Coupon Direct Yield Curve
zeroCouponDirect()
5.3.3 Zero-Coupon Bootstrapping
zeroCouponBootstrap("DataPoints")
[skipping rest of bond trading. problem is functions aren't working. strange orgbabelReoe error with Date__add__?]
6 Chapter 7: Linear Analysis
6.1 What are alpha and beta?
6.1.1 Alpha
Alpha is a measure of an investment's performance compared to a benchmark, such as the S&P 500 index.
It is an estimate of the return usually based usually on the growth of earnings per share.
6.1.2 Beta
Beta is based on the volatility, namely the movement of the fund or stock relative to its benchmark.
6.2 Linear Regression in R
idxdata <- read.csv("./inst/extdata/indices.csv", head = T, sep = ",") idxdata$Date <- as.Date(idxdata$Date, "%m/%d/%Y") head(idxdata)
Error in Date(15, "January", 2015) : could not find function "Date" Error in Date(15, "January", 2015) : could not find function "Date" Date IGSpread HYSpread SPX VIX 1 2004-09-21 53.93 363.0 1129.30 13.66 2 2004-09-22 55.18 367.5 1113.56 14.74 3 2004-09-23 54.69 371.2 1108.36 14.80 4 2004-09-24 55.11 374.5 1110.11 14.28 5 2004-09-27 54.59 370.5 1103.52 14.62 6 2004-09-28 54.15 365.5 1110.06 13.83
reg <- lm(idxdata$HYSpread ~ idxdata$SPX); reg
Call: lm(formula = idxdata$HYSpread ~ idxdata$SPX) Coefficients: (Intercept) idxdata$SPX 1542.2195 -0.7756
plot(idxdata$SPX, idxdata$HYSpread, col ="blue", main = "HY ~ SPX", xlab = "SPY", ylab = "HY") grid() abline(reg, col = "red", lw = 3)
summary(reg)
Call: lm(formula = idxdata$HYSpread ~ idxdata$SPX) Residuals: Min 1Q Median 3Q Max -361.63 -176.98 1.92 128.23 876.08 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1542.21954 23.54808 65.49 <2e-16 *** idxdata$SPX -0.77563 0.01776 -43.67 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 199.4 on 2364 degrees of freedom Multiple R-squared: 0.4465, Adjusted R-squared: 0.4462 F-statistic: 1907 on 1 and 2364 DF, p-value: < 2.2e-16
layout(matrix(1:4,2,2)) plot(reg)
reg1 <- lm(idxdata$HYSpread ~ idxdata$SPX + I(idxdata$SPX^2)) summary(reg1)
Call: lm(formula = idxdata$HYSpread ~ idxdata$SPX + I(idxdata$SPX^2)) Residuals: Min 1Q Median 3Q Max -355.6 -120.3 -34.6 130.3 422.7 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 4.399e+03 6.852e+01 64.20 <2e-16 *** idxdata$SPX -5.197e+00 1.033e-01 -50.29 <2e-16 *** I(idxdata$SPX^2) 1.658e-03 3.844e-05 43.14 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 149.1 on 2363 degrees of freedom Multiple R-squared: 0.6904, Adjusted R-squared: 0.6901 F-statistic: 2634 on 2 and 2363 DF, p-value: < 2.2e-16
plot(idxdata$SPX, idxdata$HYSpread, col = "blue", main = "HY ~ SPX", xlab = "SPX", ylab = "HY") grid() lines(idxdata$SPX, reg1$fitted.values, col = "red", lw = 3)
6.3 PCA in R
sh <- data.frame(SPX = idxdata$SPX, HY = idxdata$HYSpread) head(sh)
SPX HY 1 1129.30 363.0 2 1113.56 367.5 3 1108.36 371.2 4 1110.11 374.5 5 1103.52 370.5 6 1110.06 365.5
pca <- prcomp(sh, center = T) pca
Standard deviations (1, .., p=2): [1] 323.8234 142.0644 Rotation (n x k) = (2 x 2): PC1 PC2 SPX 0.6250655 0.7805723 HY -0.7805723 0.6250655
summary(pca)
Importance of components: PC1 PC2 Standard deviation 323.8234 142.0644 Proportion of Variance 0.8386 0.1614 Cumulative Proportion 0.8386 1.0000
source("./R/analysisHelper.R")
sh1 <- pca2dProcess(sh) head(sh1)
SPX HY PcaPrediction 1 1129.30 363.0 749.6469 2 1113.56 367.5 769.3028 3 1108.36 371.2 775.7964 4 1110.11 374.5 773.6111 5 1103.52 370.5 781.8406 6 1110.06 365.5 773.6735
plot(sh1$SPX, sh1$HY, col = "blue", main = "HY ~ SPY", xlab = "SPX", ylab = "HY") lines(sh1$SPX, sh1$PcaPrediction, type = "l", col = "red")
6.4 Comparing Linear Regression with PCA
plot(sh1$SPX, sh1$HY, col = "blue", main = "HY ~ SPY", xlab = "SPX", ylab = "HY") grid() lines(sh1$SPX, sh1$PcaPrediction, type = "l", col = "red") abline(reg, col = "darkgreen")
head(sh)
SPX HY 1 1129.30 363.0 2 1113.56 367.5 3 1108.36 371.2 4 1110.11 374.5 5 1103.52 370.5 6 1110.06 365.5
reg <- lm(sh$HY ~ sh$SPX); reg
Call: lm(formula = sh$HY ~ sh$SPX) Coefficients: (Intercept) sh$SPX 1542.2195 -0.7756
reg1 <- lm(sh$SPX ~ sh$HY); reg1
Call: lm(formula = sh$SPX ~ sh$HY) Coefficients: (Intercept) sh$HY 1610.3409 -0.5756
plot(sh$SPX, sh$HY, col = "blue", main = "HY ~ SPY", xlab = "SPX", ylab = "HY") grid() abline(reg, col = "red") lines(sh$SPX, 2797.6736 - 1.7373*sh$SPX, type = "l", col = "darkgreen")
6.5 Multiple Linear Regressions
6.6 MLR for Indices
svh <- data.frame(Date = idxdata$Date, SPX = idxdata$SPX, VIX = idxdata$VIX, HY = idxdata$HYSpread) head(svh)
Date SPX VIX HY 1 2004-09-21 1129.30 13.66 363.0 2 2004-09-22 1113.56 14.74 367.5 3 2004-09-23 1108.36 14.80 371.2 4 2004-09-24 1110.11 14.28 374.5 5 2004-09-27 1103.52 14.62 370.5 6 2004-09-28 1110.06 13.83 365.5
Call: lm(formula = svh$HY ~ svh$SPX + svh$VIX) Coefficients: (Intercept) svh$SPX svh$VIX 494.5112 -0.2735 19.3608
reg <- lm(svh$HY ~ svh$SPX + svh$VIX) reg
summary(reg)
Call: lm(formula = svh$HY ~ svh$SPX + svh$VIX) Residuals: Min 1Q Median 3Q Max -600.65 -74.31 -16.58 75.75 622.25 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 494.51124 22.22553 22.25 <2e-16 *** svh$SPX -0.27352 0.01359 -20.12 <2e-16 *** svh$VIX 19.36080 0.31123 62.21 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 122.8 on 2363 degrees of freedom Multiple R-squared: 0.7901, Adjusted R-squared: 0.79 F-statistic: 4449 on 2 and 2363 DF, p-value: < 2.2e-16
plot(svh$Date, svh$HY, type = "l", col = "blue", main = "HY ~ SPX + VIX", xlab = "Date", ylab = "HY") lines(svh$Date, reg$fitted.values, type = "l", col = "red") legend("topleft", c("HY", "Regression"), lty = c(1, 1), col = c("blue", "red")) grid()
hp <- data.frame(HY = svh$HY, MLR = reg$fitted.values) head(hp)
HY MLR 1 363.0 450.0959 2 367.5 475.3107 3 371.2 477.8947 4 374.5 467.3484 5 370.5 475.7335 6 365.5 458.6497
slr <- lm(hp$HY ~ hp$MLR) summary(slr)
Call: lm(formula = hp$HY ~ hp$MLR) Residuals: Min 1Q Median 3Q Max -600.65 -74.31 -16.58 75.75 622.25 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.992e-13 6.155e+00 0.00 1 hp$MLR 1.000e+00 1.060e-02 94.34 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 122.8 on 2364 degrees of freedom Multiple R-squared: 0.7901, Adjusted R-squared: 0.7901 F-statistic: 8901 on 1 and 2364 DF, p-value: < 2.2e-16
plot.new() plot(hp$MLR, hp$HY, col = "blue", main = "HY ~ MLR", xlab = "MLR", ylab = "HY") abline(slr, col = "red")
6.7 MLR for Stocks
library(quantmod) tks <- c("MS", "GS", "C", "XLE") stocks <- lapply(tks, function(tk) Ad(na.omit(getSymbols(tk, from = "2004-09-21", To = "2016-01-01", auto.assign = FALSE)))) stks <- do.call(merge.xts, stocks) colnames(stks) <- tks head(stks)
MS GS C XLE 2004-09-21 32.53819 77.79574 339.4682 23.75015 2004-09-22 30.26459 76.55788 330.0241 23.59989 2004-09-23 30.76155 76.29555 326.3059 23.38130 2004-09-24 31.07837 76.63162 328.3138 23.64770 2004-09-27 30.29565 76.03323 324.2982 23.63404 2004-09-28 30.34535 76.28732 329.5036 24.03705
mlr <- lm(stks$MS ~ stks$GS + stks$C + stks$XLE) mlr
Call: lm(formula = stks$MS ~ stks$GS + stks$C + stks$XLE) Coefficients: (Intercept) stks$GS stks$C stks$XLE -2.88417 0.23283 0.05075 -0.12928
reg <- lm(stks$MS ~ mlr$fitted.values) reg plot(mlr$fitted.values, stks$MS, col = "blue", main = "MS ~ MLR", xlab = "MLR", ylab = "MS") abline(reg, col = "red") grid()
Call: lm(formula = stks$MS ~ mlr$fitted.values) Coefficients: (Intercept) mlr$fitted.values 3.6e-15 1.0e+00 Error in if (length(col) < ncol(x)) col <- rep(col, length.out = ncol(x)) : argument is of length zero
6.8 Multiple PCA
isv <- data.frame(IG = idxdata$IGSpread, SPX = idxdata$SPX, VIX = idxdata$VIX) head(isv)
IG SPX VIX 1 53.93 1129.30 13.66 2 55.18 1113.56 14.74 3 54.69 1108.36 14.80 4 55.11 1110.11 14.28 5 54.59 1103.52 14.62 6 54.15 1110.06 13.83
pca <- prcomp(isv, center = T, scale. = T) pca pca$center pca$scale
Standard deviations (1, .., p=3): [1] 1.5296439 0.7327263 0.3511435 Rotation (n x k) = (3 x 3): PC1 PC2 PC3 IG 0.6008241 -0.4302517 -0.67371644 SPX -0.5087688 -0.8558801 0.09286298 VIX 0.6165750 -0.2869716 0.73313207 IG SPX VIX 89.45037 1305.46152 20.25820 IG SPX VIX 44.48189 230.79662 10.08152
pca1 <- data.frame(Pca = 0.6008241 * (isv$IG - pca$center[1])/pca$scale[1] - 0.5087688 * (isv$SPX - pca$center[2])/pca$scale[2] + 0.6165750 * (isv$VIX - pca$center[3])/pca$scale[3]) head(pca1)
Pca 1 -0.4949873 2 -0.3773544 3 -0.3688405 4 -0.3988278 5 -0.3705305 6 -0.4392060
pca1 <- data.frame(Pca=pca$x[,1]) head(pca1)
Pca 1 -0.4949873 2 -0.3773544 3 -0.3688405 4 -0.3988278 5 -0.3705305 6 -0.4392060
reg <- lm(idxdata$HYSpread ~ pca1$Pca) summary(reg)
Call: lm(formula = idxdata$HYSpread ~ pca1$Pca) Residuals: Min 1Q Median 3Q Max -387.01 -55.99 -6.30 52.22 526.28 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 529.659 1.891 280.1 <2e-16 *** pca1$Pca 164.503 1.237 133.0 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 91.99 on 2364 degrees of freedom Multiple R-squared: 0.8822, Adjusted R-squared: 0.8821 F-statistic: 1.77e+04 on 1 and 2364 DF, p-value: < 2.2e-16
plot(idxdata$Date, idxdata$HYSpread, type="l", col = "blue", main = "HY ~ IG + SPY + VIX", xlab = "Date", ylab = "HY") grid() lines(idxdata$Date, reg$fitted.values, col = "red")
reg1 <- lm(idxdata$HYSpread ~ reg$fitted.values) plot(reg$fitted.values, idxdata$HYSpread, col = "blue", main = "HY ~ IG + SPY + VIX", xlab = "PCA", ylab = "HY") abline(reg1, col = "red") grid()
6.9 PCA for Stocks
pcaData <- stks[,-1] head(pcaData)
GS C XLE 2004-09-21 77.79574 339.4682 23.75015 2004-09-22 76.55788 330.0241 23.59989 2004-09-23 76.29555 326.3059 23.38130 2004-09-24 76.63162 328.3138 23.64770 2004-09-27 76.03323 324.2982 23.63404 2004-09-28 76.28732 329.5036 24.03705
pca <- prcomp(pcaData, center = T, scale. = T) pca pca$center pca$scale
Standard deviations (1, .., p=3): [1] 1.3516952 0.9121451 0.5838763 Rotation (n x k) = (3 x 3): PC1 PC2 PC3 GS 0.5397482 0.67856504 -0.4982182 C -0.5150699 0.73432024 0.4421275 XLE 0.6658640 0.01797971 0.7458565 GS C XLE 152.38920 122.19414 51.32198 GS C XLE 45.02310 135.99251 12.28079
reg <- lm(stks$MS ~ pca$x[,1] + pca$x[,2]) summary(reg)
Call: lm(formula = stks$MS ~ pca$x[, 1] + pca$x[, 2]) Residuals: Min 1Q Median 3Q Max -9.4678 -2.2055 -0.6968 1.8374 14.3235 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 32.16287 0.05245 613.25 <2e-16 *** pca$x[, 1] 1.04637 0.03881 26.96 <2e-16 *** pca$x[, 2] 12.15232 0.05751 211.33 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 3.314 on 3989 degrees of freedom Multiple R-squared: 0.9192, Adjusted R-squared: 0.9192 F-statistic: 2.269e+04 on 2 and 3989 DF, p-value: < 2.2e-16
plot(as.Date(index(stks)), as.numeric(stks$MS), type="l", col = "blue", main = "MS ~ GS + C + XLE", xlab = "Date", ylab = "MS") grid() lines(as.Date(index(stks)), reg$fitted.values, col = "red") legend("topright", c("MS", "MLR"), lty = c(1, 1), col = c("blue", "red"))
reg1 <- lm(stks$MS ~ reg$fitted.values) plot(as.numeric(reg$fitted.values[,1]), as.numeric(stks$MS), col = "blue", main = "MS ~ GS + C + XLE", xlab = "PCA", ylab = "MS") abline(reg1, col = "red") grid()
7 Chapter 8: Time Series Analysis
7.1 Autocorrelation
\begin{equation} r_{k} = \frac{\sum_{i=1}^{N-k}(Y_{i} - \bar{Y})(Y_{i+k} - \bar{Y})} {\sum_{i=1}^{N}(Y_{i} - \bar{Y})^{2} } \end{equation}
For autocorrelation functions to make sense, the time series must be weakly stationary. Weakly stationary means that the autocorrelation is the same for any particular lag. It must satisfy these conditions:
- E(yt) is the same for all t.
- var(yt) is the same for all t.
- cov/corr between yt and yt+k is the same for all t.
7.1.1 White Noise
A collection of uncorrelated variables. If the values are drawn from a normal distribution, then its Gaussian white noise. The noise must also be iid. with mean 0 and variance finite.
set.seed(1) w<-rnorm(1000,0,1) plot.ts(w, main = "Gaussian white noise")
acf(w)
7.1.2 Random Walk with Drift
set.seed(3) w <- rnorm(300, 0, 1) y <- cumsum(w) wd <- w + 0.2 yd <- cumsum(wd) plot.ts(yd, ylim=c(-10, 80), main = "Random walk") lines(0.2*(1:300), lty="dashed") lines(y, col="red")
par(mfrow=c(2,1)) acf(y) acf(yd)
diff <- diff(y) plot.ts(diff, type="l")
acf(diff)
7.2 Model for Financial Data
library(quantmod) getSymbols("^GSPC", from = "2012-06-01", to = "2014-06-01") plot(Cl(GSPC))
acf(Cl(GSPC))
rtn <- diff(log(Cl(GSPC))) plot(rtn)
acf(rtn,na.action=na.pass)
7.2.1 MOdel Selection
- AIC (Akaike Information Criterion): 2k - 2ln(
- AICc takes into account sample size: 2kn/(n-k-1)-2ln(L)
- BIC: -2ln(L) + kln(n)
AIC assumes true model infinite number of parameters and attempts to minimize information loss using finite dimensional model. As n increases, selection based on AIC equals that based on maximum likelihood. AIC favors more complex models with more data. BIC is a large-sample approximation to Bayesian selection among fixed set of finite dimensional models. BIC does not become equivalent to maximum likelihood selection. BIC will favor consistent selection among a fixed family of models.
7.2.2 ARIMA Model
- AR Model (autoregressive model)
Data needs to be stationary. Timepoint y is predicted by alpha * previous timepoint + w (error). Stationary is when all absolute values of the roots of equation is greater than unity. The equation is the characteristic equation: 1 - alpha1*B - alpha2*B2 - … - alphap*Bp = 0 Random walk is an AR(1) model when alpha = 1.
7.2.3 MA model (moving average model)
linear combination of past white noise terms (with some beta) Means are 0.
set.seed(2) y <- arima.sim(n = 100, list(ma = 0.5)) layout(1:2) plot(y) acf(y)
auto.arima(y)
Series: y ARIMA(0,0,1) with zero mean Coefficients: ma1 0.4219 s.e. 0.0867 sigma^2 estimated as 1.346: log likelihood=-156.35 AIC=316.7 AICc=316.83 BIC=321.92
set.seed(1) y <- arima.sim(n = 200, list(ma = c(0.7, 0.5))) layout(1:2) plot(y) ic <- acf(y)
auto.arima(y,max.p=0)
Series: y ARIMA(0,0,2) with zero mean Coefficients: ma1 ma2 0.6735 0.5145 s.e. 0.0625 0.0596 sigma^2 estimated as 0.8782: log likelihood=-270.21 AIC=546.43 AICc=546.55 BIC=556.32
7.2.4 Autoregressive Moving Average Model (ARMA)
Past values and white noise. Special form of ARIMA, without differencing.
set.seed(1) y <- arima.sim(n = 200, list(ar = 0.6, ma = 0.3)) layout(1:2) plot(y) acf(y)
arima(y,order=c(1,0,1))
Call: arima(x = y, order = c(1, 0, 1)) Coefficients: ar1 ma1 intercept 0.6094 0.2587 0.1367 s.e. 0.0745 0.0913 0.2130 sigma^2 estimated as 0.8882: log likelihood = -272.34, aic = 552.68
set.seed(1) y <- arima.sim(n = 500, list(ar = c(0.8, -0.5), ma = c(0.6, -0.4))) layout(1:2) plot(y, type = "l") acf(y)
arima(y,order=c(2,0,2))
Call: arima(x = y, order = c(2, 0, 2)) Coefficients: ar1 ar2 ma1 ma2 intercept 0.8283 -0.5139 0.5488 -0.4512 0.0060 s.e. 0.0761 0.0392 0.0871 0.0869 0.0726 sigma^2 estimated as 1.027: log likelihood = -719.95, aic = 1451.91
set.seed(2) y <- arima.sim(n=1000, list(ar = c(0.3, 0.9, 0.1, -0.8), ma = c(0.9, 0.8, 0.72))) layout(1:2) plot(y, type = "l") acf(y)
fit<-auto.arima(y) fit
Series: y ARIMA(4,0,3) with non-zero mean Coefficients: ar1 ar2 ar3 ar4 ma1 ma2 ma3 mean 0.3314 0.8765 0.0906 -0.7818 0.8996 0.8270 0.7257 0.4081 s.e. 0.0208 0.0244 0.0241 0.0206 0.0243 0.0236 0.0248 0.2211 sigma^2 estimated as 0.9656: log likelihood=-1403.22 AIC=2824.44 AICc=2824.62 BIC=2868.61
acf(resid(fit))
The Ljung-Box test confirms the white noise model of the residuals.
Box.test(resid(fit),lag=30,type="Ljung-Box")
getSymbols("^GSPC", from = "1980-01-01", to = "2016-01-01") spx <- na.omit(diff(log(Ad(GSPC)))) layout(1:2) plot(spx, type="l") acf(spx)
fit <- auto.arima(spx) fit Box.test(resid(fit), lag=30, type = "Ljung-Box")
Series: spx ARIMA(2,0,1) with non-zero mean Coefficients: ar1 ar2 ma1 mean 0.6523 -0.0188 -0.6807 3e-04 s.e. 0.1118 0.0133 0.1114 1e-04 sigma^2 estimated as 0.0001263: log likelihood=27871.21 AIC=-55732.43 AICc=-55732.42 BIC=-55696.86 Box-Ljung test data: resid(fit) X-squared = 77.632, df = 30, p-value = 4.302e-06
There is poor fit. In practice, ARMA models are not good models for financial data.
7.2.5 ARIMA Model Applications
set.seed(1) y <- arima.sim(n = 1000, list(order = c(1,1,1), ar = 0.8, ma = -0.6)) layout(1:2) plot(y) acf(y)
Use auto arima and AIC criterion.
aic.fit <- auto.arima(y, approximation = FALSE, ic = "aic") aic.fit
Series: y ARIMA(2,1,2) Coefficients: ar1 ar2 ma1 ma2 -0.1073 0.7623 0.2801 -0.6682 s.e. 0.0724 0.0604 0.0857 0.0793 sigma^2 estimated as 1.064: log likelihood=-1448.28 AIC=2906.55 AICc=2906.61 BIC=2931.09
Using BIC, we recover ARIMA(1,1,1)
bic.fit <- auto.arima(y, approximation = FALSE, ic = "bic") bic.fit
Series: y ARIMA(1,1,1) Coefficients: ar1 ma1 0.7972 -0.6502 s.e. 0.0714 0.0912 sigma^2 estimated as 1.075: log likelihood=-1454.07 AIC=2914.13 AICc=2914.16 BIC=2928.86
Run time series diagnostics.
tsdiag(aic.fit)
tsdiag(bic.fit)
Checking fit
plot(y) lines(fitted(aic.fit), col = "red") lines(fitted(bic.fit), col = "green")
ARIMA applied to SPX
library(quantmod) getSymbols("^GSPC", from = "2007-01-01", to = "2016-01-01") sp <- Ad(GSPC) plot(sp)
sp.arima <- auto.arima(sp, approximation = FALSE) sp.arima
Series: sp ARIMA(1,1,1) Coefficients: ar1 ma1 0.5726 -0.6529 s.e. 0.1460 0.1354 sigma^2 estimated as 272.1: log likelihood=-9561.89 AIC=19129.78 AICc=19129.79 BIC=19146.96
Residual plot
plot(sp.arima$residuals)
Other diagnostics with tsdiag
tsdiag(sp.arima)
Forecast with sp arima model
plot(forecast(sp.arima, h = 30))
7.3 GARCH Model
In ARIMA models, conditional variance (past variance) is constant. However, there are clusters of variability in the SPX residuals, suggesting variance is not constant. ARIMA models provide linear forecasts. Need to use generalized autoregressive conditional heteroscdasticity (GARCH) models.
To see if GARCH is needed, check for clusters of variability in residuals and squared residuals. ACF of squared residuals also help confirm whether residuals are dependent.
Let's check squared residuals acf.
res2 <- sp.arima$residuals^2 layout(1:2) plot(res2) acf(res2)
ARCH assumes that error variance (sometimes called innovation) is related to previous squared error. It's like applying AR(q) model to error variance.
Financial time series are usually modeled in two phases. First, find best fitting ARIMA. Plot residuals and look for nonlinearity or variabliity clusters. Then use GARCH to model the nonlinear patterns of the residuals. Then apply the hybrid model, combinging ARIMA and GARCH.
sp.garch <- garch(sp.arima$residuals, order = c(1, 1), trace = FALSE) summary(sp.garch)
Call: garch(x = sp.arima$residuals, order = c(1, 1), trace = FALSE) Model: GARCH(1,1) Residuals: Min 1Q Median 3Q Max -5.6793 -0.4959 0.1168 0.6432 3.0569 Coefficient(s): Estimate Std. Error t value Pr(>|t|) a0 7.32201 1.24338 5.889 3.89e-09 *** a1 0.09988 0.01025 9.748 < 2e-16 *** b1 0.87142 0.01363 63.938 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Diagnostic Tests: Jarque Bera Test data: Residuals X-squared = 262.88, df = 2, p-value < 2.2e-16 Box-Ljung test data: Squared.Residuals X-squared = 3.9171, df = 1, p-value = 0.0478
Test goodness of fit by plotting GARCH residuals.
layout(1:2) acf(sp.garch$residuals, na.action = na.pass) acf(sp.garch$residuals^2, na.action = na.pass)
The next step would be to get the forecasts from this, but the book doesn't show how to do that.
7.4 Cointegration
The idea is the find a linear combination of two or more time series that yields a stationary time series. One finance strategy is to use a mean-reerting strategy, where buy low, wait for reversion to mean, then sell at the higher price. However, most price series are not mean-reverting. We can apply a mean-reverting strategy if we make the time series stationary.
7.4.1 Testing for Stationarity
We can test for stationarity using the Augmented Dickey Fuller (ADF) test. The ADF testswhether capital lambda is 0, which indicates that price changes are completely independent.
library(quantmod) getSymbols("SPY", to = "2016-01-01") spy <- Ad(SPY) plot.new() plot(spy)
Augmented Dickey Fuller Test for SPY
spy.adf <- ur.df(spy, type = "trend", selectlags = "AIC") summary(spy.adf)
############################################### # Augmented Dickey-Fuller Test Unit Root Test # ############################################### Test regression trend Call: lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag) Residuals: Min 1Q Median 3Q Max -7.9014 -0.6377 0.0781 0.7021 9.9537 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.566e-01 1.139e-01 1.375 0.1693 z.lag.1 -2.643e-03 1.377e-03 -1.919 0.0551 . tt 1.796e-04 7.564e-05 2.374 0.0177 * z.diff.lag -5.149e-02 2.101e-02 -2.451 0.0143 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.358 on 2260 degrees of freedom Multiple R-squared: 0.005144, Adjusted R-squared: 0.003824 F-statistic: 3.895 on 3 and 2260 DF, p-value: 0.008649 Value of test-statistic is: -1.9187 2.4312 2.8176 Critical values for test statistics: 1pct 5pct 10pct tau3 -3.96 -3.41 -3.12 phi2 6.09 4.68 4.03 phi3 8.27 6.25 5.34
The book says that tau3 = -1.9187, phi2 = 2.4312, and phi3=2.9176. The book says we fail to reject the null for phi3, indicating SPY daily price contains a unit root. This is confirmed by tar3, which is not more negative than 10% level. So we cannot reject a unit root. Then we examine whether SPY is a random walk with or without drift–if so, then not stationary. phi2 is less than critical values, so SPY acts as a random walk. there are more test. The book suggests we read more about them. (p.285)
7.4.2 Cointegrating time Series
There are two approaches.
- covariate Augmented Dickery-Fuller Test
First approach is the covariate (cointegration) Augmented Dickey-Fuller Test. First proposed by Eagle and Granger (1987), tis approach performs lr to fit the two or more time series. Then you apply the normal ADF test.
SPY and IWM Example
library(quantmod) getSymbols("IWM", to = "2016-01-01") iwm <- Ad(IWM) plot.new() plot(spy, ylim = c(30,210), main = ("SPY and IWM")) lines(iwm, col = "red") #legend("topleft", c("SPY", "IWM"), col = c("black", "red"), lty = c(1,1))
Regress spy on iwm
lm(spy ~ iwm + 1)
Call: lm(formula = spy ~ iwm + 1) Coefficients: (Intercept) iwm 4.111 1.586
Test for stationarity of residual series:
res <- ur.df(spy-4.468-1.626*iwm, type = "none", selectlags = "AIC") summary(res)
############################################### # Augmented Dickey-Fuller Test Unit Root Test # ############################################### Test regression none Call: lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag) Residuals: Min 1Q Median 3Q Max -5.0930 -0.4494 -0.0340 0.4162 3.7131 Coefficients: Estimate Std. Error t value Pr(>|t|) z.lag.1 -0.005789 0.002391 -2.422 0.0155 * z.diff.lag -0.025443 0.021030 -1.210 0.2265 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.7266 on 2262 degrees of freedom Multiple R-squared: 0.003399, Adjusted R-squared: 0.002518 F-statistic: 3.857 on 2 and 2262 DF, p-value: 0.02127 Value of test-statistic is: -2.4217 Critical values for test statistics: 1pct 5pct 10pct tau1 -2.58 -1.95 -1.62
-2.4217 is less than 5%, so we can reject null that lambda is zero. That is, SPY and IWM are cointegrating at the .05 level.
- Johansen Test
Our goal is the form a stationary price series out of an arbirary number of series. We then examine whether capital lambda (a matrix) can be rejected. r equals the number of independently combined series. Johansen computer r baed on eigenvector decomposition of capital lambda to test whether we can reject the null hypothesis that r = n. The eigenvectors are hedge ratios and can then be combined to form a stationary series.
Example
Looks like we mean normalize the scores first.
getSymbols(c("MS", "GS", "C"), from = "2010-01-01", to ="2016-01-01") ms <- Ad(MS) gs <- Ad(GS) c <- Ad(C) plot(ms/mean(ms), main = "Mean normalized prices for MS, GS, and C") lines(gs/mean(gs), col = "red") lines(c/mean(c), col = "blue") #legend("topleft", c("MS", "GS", "C"), col = c("black", "red", "blue"), lty = c(1,1,1))
Perform Johansen test on three stocks
prices <- na.omit(merge(ms, gs, c)) result <- ca.jo(prices, type = "eigen", ecdet = "const") summary(result)
###################### # Johansen-Procedure # ###################### Test type: maximal eigenvalue statistic (lambda max) , without linear trend and constant in cointegration Eigenvalues (lambda): [1] 2.642826e-02 8.175185e-03 9.652049e-04 2.658500e-19 Values of teststatistic and critical values of test: test 10pct 5pct 1pct r <= 2 | 1.46 7.52 9.24 12.97 r <= 1 | 12.38 13.75 15.67 20.20 r = 0 | 40.39 19.77 22.00 26.81 Eigenvectors, normalised to first column: (These are the cointegration relations) MS.Adjusted.l2 GS.Adjusted.l2 C.Adjusted.l2 constant MS.Adjusted.l2 1.0000000 1.0000000 1.0000000 1.000000 GS.Adjusted.l2 -0.4338740 -0.1573598 0.3349882 1.251625 C.Adjusted.l2 0.8562129 -0.2441210 -0.2506577 -19.160619 constant 2.7190780 8.3617592 -62.1702500 115.074221 Weights W: (This is the loading matrix) MS.Adjusted.l2 GS.Adjusted.l2 C.Adjusted.l2 constant MS.Adjusted.d -0.001238315 -0.012044880 -0.0008156064 2.020359e-19 GS.Adjusted.d 0.021507309 0.008929769 -0.0044987255 5.237316e-19 C.Adjusted.d -0.026951439 -0.002360268 -0.0012925429 3.771185e-19
We can only reject null when r = 0. We construct the stationary series using the first column, which corresponds to the largest eigenvalue. (is this true?)
Plot combined stationary series
spread = ms - 0.4517*gs + 0.8659*c + 2.8859 plot(spread)
It looks like the combined series showed mean-reverting behavior.
8 Chapter 9: Machine Learning
Supervised learning is widely used in finance.
8.1 KNN Classifier
Non-parametric and can be used for classification or regression.
8.1.1 KNN Model
Near objects have similar characteristics. New instances classified by majority vote of k-neighbors. Based on Euclidiean distance.
- KNN Model for Classification
In this section, we will predict stock price direction. We will use global stock data as input features to predict moving direction for FTSE 100 index. We believe that the timing of Nikkei 225 (Tokyo) influences London stock prices the next day. IVs: open, high, low, close, and daily returns for FTSE and N225. DV: today's market move direction, up or down (+ or -_
We can use todays N225 as predictors because Tokyo sotck market already closed before London market opens. However, we cannot use today's direction for London because today hasn't happened yet.
Yahoo does not provie recent historical data for FTSE. Se we load data from a csv file. Use the lead command to offset the direction field by one day. (yesterday -> today)
Load Machinelearning.R file Load FTSE and N225 data
fn <- processFtseN225Data() head(fn, 10)
Date FTSE.Open FTSE.High FTSE.Low FTSE.Close FTSE.Rtn N225.Open 1 2010-01-04 5412.88 5500.34 5410.82 5500.34 0.00000000 10719.44 2 2010-01-05 5500.34 5536.38 5480.71 5522.50 0.40288419 10709.55 3 2010-01-06 5522.50 5536.48 5497.65 5530.04 0.13653237 10742.75 4 2010-01-07 5530.04 5551.66 5499.80 5526.72 -0.06003573 10743.30 5 2010-01-08 5526.72 5549.25 5494.79 5534.24 0.13606624 10770.35 6 2010-01-12 5538.07 5549.61 5459.94 5498.71 -0.71071691 10795.48 7 2010-01-13 5498.71 5509.67 5450.85 5473.48 -0.45883489 10778.07 8 2010-01-14 5473.48 5521.90 5473.48 5498.20 0.45163223 10917.41 9 2010-01-15 5498.20 5527.02 5450.38 5455.37 -0.77898221 10887.61 10 2010-01-18 5455.37 5504.00 5454.31 5494.39 0.71525854 10866.83 N225.High N225.Low N225.Close N225.Rtn Direction 1 10791.04 10655.57 10681.83 0.2537829 UP 2 10768.61 10661.17 10731.45 0.4645282 UP 3 10774.00 10636.67 10681.66 -0.4639638 DOWN 4 10816.45 10677.56 10798.32 1.0921538 UP 5 10905.39 10763.68 10879.14 0.7484436 DOWN 6 10866.62 10729.86 10735.03 -1.3246394 DOWN 7 10909.94 10774.25 10907.68 1.6082807 UP 8 10982.10 10878.83 10982.10 0.6822709 DOWN 9 10895.10 10781.03 10855.08 -1.1566052 UP 10 10866.83 10749.47 10764.90 -0.8307602 UP
Split data into test and train
fn.train <- with(fn, fn[Date >= "2010-01-01" & Date <= "2014-12-31",]) fn.test <- with(fn, fn[Date >= "2015-01-01",])
Train the KNN Model using the train function
set.seed(1) fn.knn <- train(Direction ~ ., data = fn.train[,-1], method = "knn", trControl = trainControl(method = "repeatedcv", repeats = 3), preProcess = c("center","scale"), tuneLength = 50)
This model uses 10-fold crossvalidation, repeated (here) 3 times. It would be good to learn how it does this from the documentation.
Show model results
fn.knn
k-Nearest Neighbors 1192 samples 10 predictor 2 classes: 'DOWN', 'UP' Pre-processing: centered (10), scaled (10) Resampling: Cross-Validated (10 fold, repeated 3 times) Summary of sample sizes: 1073, 1072, 1072, 1073, 1073, 1073, ... Resampling results across tuning parameters: k Accuracy Kappa 5 0.5406026 0.07707570 7 0.5307705 0.05757382 9 0.5380489 0.07159595 11 0.5495243 0.09330675 13 0.5525892 0.09863465 15 0.5512003 0.09538997 17 0.5660395 0.12559369 19 0.5705001 0.13431845 21 0.5651873 0.12368487 23 0.5663005 0.12511593 25 0.5704837 0.13365812 27 0.5796830 0.15231048 29 0.5867000 0.16637323 31 0.5881076 0.16921754 33 0.5889434 0.17101132 35 0.5920388 0.17754017 37 0.5951129 0.18357861 39 0.5973421 0.18759681 41 0.5931404 0.17873537 43 0.5981754 0.18911737 45 0.6006848 0.19395187 47 0.6028953 0.19854030 49 0.6015156 0.19550108 51 0.5992654 0.19077423 53 0.6031799 0.19839196 55 0.5998232 0.19138709 57 0.6006660 0.19307620 59 0.6051456 0.20210594 61 0.6009462 0.19358572 63 0.6026268 0.19606696 65 0.6031752 0.19691762 67 0.6040133 0.19820572 69 0.6073677 0.20457490 71 0.6026196 0.19486267 73 0.6023395 0.19455030 75 0.6023301 0.19448917 77 0.5973093 0.18420313 79 0.5989806 0.18767030 81 0.5961795 0.18204971 83 0.5956404 0.18082816 85 0.5989924 0.18712784 87 0.5953461 0.17981105 89 0.5945198 0.17811865 91 0.5922788 0.17344317 93 0.5922764 0.17333611 95 0.5939243 0.17678026 97 0.5939313 0.17685368 99 0.5875169 0.16408313 101 0.5877876 0.16476684 103 0.5889034 0.16700159 Accuracy was used to select the optimal model using the largest value. The final value used for the model was k = 69.
Show accuracy
plot(fn.knn)
Predict on test set with confusion matrix
fn.knn.predict <- predict(fn.knn, newdata = fn.test[,-1]) confusionMatrix(fn.knn.predict, as.factor(fn.test$Direction))
Confusion Matrix and Statistics Reference Prediction DOWN UP DOWN 100 76 UP 16 45 Accuracy : 0.6118 95% CI : (0.5466, 0.6742) No Information Rate : 0.5105 P-Value NIR] : 0.001082 Kappa : 0.2315 Mcnemar's Test P-Value : 7.691e-10 Sensitivity : 0.8621 Specificity : 0.3719 Pos Pred Value : 0.5682 Neg Pred Value : 0.7377 Prevalence : 0.4895 Detection Rate : 0.4219 Detection Prevalence : 0.7426 Balanced Accuracy : 0.6170 'Positive' Class : DOWN
fn.test$Predict <- fn.knn.predict head(fn.test)
Date FTSE.Open FTSE.High FTSE.Low FTSE.Close FTSE.Rtn N225.Open 1193 2015-01-05 6547.80 6576.74 6404.49 6417.16 -1.9951740 17101.58 1194 2015-01-06 6417.16 6452.66 6328.59 6366.51 -0.7892900 16808.26 1195 2015-01-07 6366.51 6459.74 6366.51 6419.83 0.8375075 17067.40 1196 2015-01-08 6419.83 6580.82 6419.83 6569.96 2.3385354 17318.74 1197 2015-01-09 6569.96 6570.24 6471.38 6501.14 -1.0474950 16970.88 1198 2015-01-13 6501.42 6558.83 6465.19 6542.20 0.6272476 16961.82 N225.High N225.Low N225.Close N225.Rtn Direction Predict 1193 17111.36 16881.73 16883.19 -3.01872716 DOWN DOWN 1194 16974.61 16808.26 16885.33 0.01267903 UP DOWN 1195 17243.71 17016.09 17167.10 1.66872385 UP UP 1196 17342.65 17129.53 17197.73 0.17842769 DOWN UP 1197 17087.71 16828.27 17087.71 -0.63973285 UP DOWN 1198 17036.72 16770.56 16795.96 -1.70736737 DOWN DOWN
Compare predicted results with real market moves
fn.test$Predict <- fn.knn.predict head(fn.test, 20)
8.1.2 Confusion Matrix
Accuracy is along diagonal. Easy to calculate.
8.1.3 KNN Model for Regression
Change Direction to Expected to predict continuous values
fnReg <- fn colnames(fnReg)[colnames(fnReg)=="Direction"] <- "Expected"
Offset the FTSE close prices by one day and place them into the Expected column
fnReg$Expected <- lead(fnReg$FTSE.Close) fnReg <- fnReg[-nrow(fnReg),] #Remove the last row with NA head(fnReg)
Date FTSE.Open FTSE.High FTSE.Low FTSE.Close FTSE.Rtn N225.Open 1 2010-01-04 5412.88 5500.34 5410.82 5500.34 0.00000000 10719.44 2 2010-01-05 5500.34 5536.38 5480.71 5522.50 0.40288419 10709.55 3 2010-01-06 5522.50 5536.48 5497.65 5530.04 0.13653237 10742.75 4 2010-01-07 5530.04 5551.66 5499.80 5526.72 -0.06003573 10743.30 5 2010-01-08 5526.72 5549.25 5494.79 5534.24 0.13606624 10770.35 6 2010-01-12 5538.07 5549.61 5459.94 5498.71 -0.71071691 10795.48 N225.High N225.Low N225.Close N225.Rtn Expected 1 10791.04 10655.57 10681.83 0.2537829 5522.50 2 10768.61 10661.17 10731.45 0.4645282 5530.04 3 10774.00 10636.67 10681.66 -0.4639638 5526.72 4 10816.45 10677.56 10798.32 1.0921538 5534.24 5 10905.39 10763.68 10879.14 0.7484436 5498.71 6 10866.62 10729.86 10735.03 -1.3246394 5473.48
Use the last 60 days to split into train and test.
fnReg.train <- head(fnReg, -60) fnReg.test <- tail(fnReg, 60)
Train the Knn model
set.seed(1) fnReg.knn <- train(Expected ~ ., data = fnReg.train[, -1], method = "knn", trControl = trainControl(method = "repeatedcv", repeats = 3), tuneLength = 20)
Look at the output
fnReg.knn
Error: object 'fnReg.knn' not found
Plot results to see the fit [this plot isn't working]
plot(fnReg.train$Expected, type="l", ylab = "FTSE Price") lines(predict(fnReg.knn, newdata = fnReg.train[, -1]), col="red") #legend("topleft", c("Training Data", "Prediction"), col = c("black", "red"), lty = c(1, 1))
Use trained model to predict data
fnReg.predict <- predict(fnReg.knn, newdata = fnReg.test[, -1])
Error in predict(fnReg.knn, newdata = fnReg.test[, -1]) : object 'fnReg.knn' not found
Use default summary to evaluate performance of KNN model.
modelValues <- data.frame(obs = fnReg.test$Expected, pred = fnReg.predict) defaultSummary(modelValues)
Error in data.frame(obs = fnReg.test$Expected, pred = fnReg.predict) : object 'fnReg.predict' not found Error in defaultSummary(modelValues) : object 'modelValues' not found
Plot test data and prediction
plot(fnReg.test$Expected, type = "l", ylim = c(4500, 7500), ylab = "FTSE Price") lines(fnReg.predict, col = "red") #legend("topright", c("Test Data", "Prediction"), col = c("black", "red"), #lty = c(1, 1))
8.1.4 Random Forest Algorithm
Ensemble weak learners to make a strong learner.
- Take a random sample of size N with replacement from the data (bootstrap sample).
- Take a random sample without replacements of the predictors.
- Construct a split by using predictors selected in Step 2.
- Repeat Steps 2 and 3 for each subsequent split until the tree is as large as desired.
- Drop the out-of-bag data down the tree. Store the class assigned to each observation along with each observation's predictor value.
- Repeat Steps 1-5 a large number of time (e.g., 1000).
- For each observation in the dataset, count the number of trees that it is classified in one category over the number of trees.
- Assign each obsevation to a final category by a majority vote over the set of trees.
Bootstrapping decreases variance without increasing bias. Optimal number of trees can be found using cross-validation.
8.1.5 Random Forest for Classification
Divide into test and train
fn.train <- with(fn, fn[Date <= "2014-12-31",]) fn.test <- with(fn, fn[Date >= "2015-01-01",])
Train using random forest
set.seed(1) fn.rf <- train(Direction ~ ., data = fn.train[,-1], method = "rf", trControl = trainControl(method = "repeatedcv", repeats = 3), preProcess = c("center","scale"), tuneLength = 50) fn.rf
note: only 9 unique complexity parameters in default grid. Truncating the grid to 9 . Random Forest 1192 samples 10 predictor 2 classes: 'DOWN', 'UP' Pre-processing: centered (10), scaled (10) Resampling: Cross-Validated (10 fold, repeated 3 times) Summary of sample sizes: 1073, 1072, 1072, 1073, 1073, 1073, ... Resampling results across tuning parameters: mtry Accuracy Kappa 2 0.5654930 0.1254466 3 0.5680329 0.1306635 4 0.5716605 0.1381195 5 0.5713921 0.1378262 6 0.5736309 0.1421767 7 0.5744499 0.1432106 8 0.5800572 0.1552218 9 0.5783766 0.1515717 10 0.5802974 0.1556348 Accuracy was used to select the optimal model using the largest value. The final value used for the model was mtry = 10.
Plot accuracy y
plot(fn.rf)
Predict the test set
fn.rf.predict <- predict(fn.rf, newdata = fn.test[,-1]) confusionMatrix(fn.rf.predict,as.factor(fn.test$Direction))
Error in predict(fn.rf, newdata = fn.test[, -1]) : object 'fn.rf' not found Error in confusionMatrix(fn.rf.predict, as.factor(fn.test$Direction)) : object 'fn.rf.predict' not found
8.1.6 Random Forest for Regression
Train on continuous variables
set.seed(1) fnReg.rf <- train(Expected ~ ., data = fnReg.train[, -1], method = "rf", trControl = trainControl(method = "repeatedcv", repeats = 3), tuneLength = 20)
note: only 9 unique complexity parameters in default grid. Truncating the grid to 9 .
Preview results
fnReg.rf
Random Forest 1368 samples 10 predictor No pre-processing Resampling: Cross-Validated (10 fold, repeated 3 times) Summary of sample sizes: 1232, 1230, 1232, 1230, 1232, 1232, ... Resampling results across tuning parameters: mtry RMSE Rsquared MAE 2 58.85740 0.9888640 44.01967 3 58.65095 0.9889508 43.82729 4 58.77135 0.9889082 43.94228 5 58.94953 0.9888421 44.07535 6 59.14551 0.9887648 44.22918 7 59.33560 0.9886933 44.32935 8 59.47804 0.9886415 44.49377 9 59.55471 0.9886139 44.57188 10 59.68082 0.9885666 44.67152 RMSE was used to select the optimal model using the smallest value. The final value used for the model was mtry = 3.
Predict random forest regression
fnReg.predict <- predict(fnReg.rf, newdata = fnReg.test[, -1]) modelValues <- data.frame(obs = fnReg.test$Expected, pred = fnReg.predict) defaultSummary(modelValues)
RMSE Rsquared MAE 97.6879629 0.5938378 82.7876686
Plot results
plot(fnReg.test$Expected, type = "l", ylim = c(4500, 7500), ylab = "FTSE Price") lines(fnReg.predict, col = "red") legend("topright", c("Test Data", "Prediction"), col = c("black", "red"), lty = c(1, 1))
8.2 Support Vector Machine
Finds the hypoerplan that gives the largest minimum distance to the data points in the training data set. This distance is called margin. Therefore, SVM maximizes margin. SVM also works for nonlinear classification.
The large-margin separations means that only the relative position or similarlity of the data points to each other is important. Similarity can be computered using a kernel function. The simplest kernel function is the dot-product between two feature vectors (linear kernel). Polynomial and Gaussian kernels are also possible, and these are nonlinear.
The loss function in SVM should be modified to include a distance measure, which is often called epsilon intensive loss function. This epsilon intensive loss functions ensures existence of a global minimum.
There are two steps. First, transform input data to high-dimensional feature space by specifying a kernel function. The step is to solve a quadratic optimization problem to fit optimal hypoterplane to classify transformed features into different classes. The number of transformed features is determined by the number of support vectors.
8.2.1 SVM for Classification
Predict next day's stock movement.
ctrl <- trainControl(method = "repeatedcv", repeats = 10) fn.svm <- train(Direction ~ ., data = fn.train[,-1], method = "svmLinear", trControl = ctrl, preProcess = c("center","scale")) fn.svm fn.svm.predict <- predict(fn.svm, newdata = fn.test[,-1]) confusionMatrix(fn.svm.predict,as.factor(fn.test$Direction))
Support Vector Machines with Linear Kernel 1192 samples 10 predictor 2 classes: 'DOWN', 'UP' Pre-processing: centered (10), scaled (10) Resampling: Cross-Validated (10 fold, repeated 10 times) Summary of sample sizes: 1072, 1072, 1073, 1073, 1073, 1073, ... Resampling results: Accuracy Kappa 0.5995895 0.1927101 Tuning parameter 'C' was held constant at a value of 1 Confusion Matrix and Statistics Reference Prediction DOWN UP DOWN 34 16 UP 82 105 Accuracy : 0.5865 95% CI : (0.5209, 0.6499) No Information Rate : 0.5105 P-Value NIR] : 0.01132 Kappa : 0.1628 Mcnemar's Test P-Value : 5.169e-11 Sensitivity : 0.2931 Specificity : 0.8678 Pos Pred Value : 0.6800 Neg Pred Value : 0.5615 Prevalence : 0.4895 Detection Rate : 0.1435 Detection Prevalence : 0.2110 Balanced Accuracy : 0.5804 'Positive' Class : DOWN
8.2.2 SVM for Regression
Train, view results, predict test
set.seed(1234) fnReg.svm <- train(Expected ~ ., data = fnReg.train[, -1], method = "svmLinear", trControl = trainControl(method = "repeatedcv", repeats = 3)) fnReg.svm fnReg.predict <- predict(fnReg.svm, newdata = fnReg.test[, -1]) modelValues <- data.frame(obs = fnReg.test$Expected, pred = fnReg.predict) defaultSummary(modelValues)
Support Vector Machines with Linear Kernel 1368 samples 10 predictor No pre-processing Resampling: Cross-Validated (10 fold, repeated 3 times) Summary of sample sizes: 1232, 1230, 1230, 1232, 1230, 1232, ... Resampling results: RMSE Rsquared MAE 57.15723 0.9895976 42.67837 Tuning parameter 'C' was held constant at a value of 1 RMSE Rsquared MAE 64.1340012 0.8099375 49.8210808
Plot results
plot(fnReg.test$Expected, type = "l", ylim = c(4500, 7500), ylab = "FTSE Price") lines(fnReg.predict, col = "red") #legend("topright", c("Test Data", "Prediction"), col = c("black", "red"), #lty = c(1, 1))
8.3 Artificial Neural Networks (ANNs)
This book suggests that multilayer feedforward ANNs are best for predicting stock price movemement (no connections backward).
8.3.1 Neural Networks for Classification
First we will examine the data to see if we should normalize our features (the answer is always yes)
str(fn)
'data.frame': 1429 obs. of 12 variables: $ Date : Date, format: "2010-01-04" "2010-01-05" ... $ FTSE.Open : num 5413 5500 5522 5530 5527 ... $ FTSE.High : num 5500 5536 5536 5552 5549 ... $ FTSE.Low : num 5411 5481 5498 5500 5495 ... $ FTSE.Close: num 5500 5522 5530 5527 5534 ... $ FTSE.Rtn : num 0 0.403 0.137 -0.06 0.136 ... $ N225.Open : num 10719 10710 10743 10743 10770 ... $ N225.High : num 10791 10769 10774 10816 10905 ... $ N225.Low : num 10656 10661 10637 10678 10764 ... $ N225.Close: num 10682 10731 10682 10798 10879 ... $ N225.Rtn : num 0.254 0.465 -0.464 1.092 0.748 ... $ Direction : chr "UP" "UP" "DOWN" "UP" ...
Now we will create a function to normalize the data to the range.
normalizeData <- function(x) { return((x-min(x)) / (max(x)-min(x))) } fnx <- fn fnx$Date <- fnx$Direction <- NULL fnx <- as.data.frame(lapply(fnx, normalizeData)) head(fnx)
FTSE.Open FTSE.High FTSE.Low FTSE.Close FTSE.Rtn N225.Open N225.High 1 0.2641729 0.2819301 0.2736763 0.3022282 0.4748689 0.2019377 0.2031963 2 0.3022282 0.2978813 0.3044879 0.3118704 0.5158622 0.2011639 0.2014375 3 0.3118704 0.2979255 0.3119561 0.3151512 0.4887610 0.2037613 0.2018601 4 0.3151512 0.3046442 0.3129039 0.3137066 0.4687603 0.2038043 0.2051887 5 0.3137066 0.3035775 0.3106952 0.3169787 0.4887136 0.2059205 0.2121627 6 0.3186452 0.3037368 0.2953313 0.3015190 0.4025538 0.2078865 0.2091227 N225.Low N225.Close N225.Rtn 1 0.1983699 0.1984432 0.5917886 2 0.1988108 0.2023478 0.6033282 3 0.1968820 0.1984298 0.5524876 4 0.2001010 0.2076099 0.6376946 5 0.2068808 0.2139696 0.6188743 6 0.2042184 0.2026296 0.5053603
Split data into test and train
x.train <- head(fnx, -241) x.test <- tail(fnx, 241) y.train <- as.data.frame(head(fn,-241))$Direction y.test <- as.data.frame(tail(fn,241))$Direction
Train ANN
set.seed(1234) ctrl <- trainControl(method = "repeatedcv", repeats = 5, number = 5) fn.nn <- train(x = x.train, y = y.train, method = "nnet", trControl = ctrl, tuneLength = 10, trace = FALSE) fn.nn
Neural Network 1188 samples 10 predictor 2 classes: 'DOWN', 'UP' No pre-processing Resampling: Cross-Validated (5 fold, repeated 5 times) Summary of sample sizes: 950, 951, 950, 951, 950, 950, ... Resampling results across tuning parameters: size decay Accuracy Kappa 1 0.0000000000 0.5947506 0.1773183 1 0.0001000000 0.5840008 0.1489377 1 0.0002371374 0.5767667 0.1340513 1 0.0005623413 0.5645033 0.1063615 1 0.0013335214 0.5924027 0.1721731 1 0.0031622777 0.5851829 0.1546109 1 0.0074989421 0.5967823 0.1837368 1 0.0177827941 0.5961072 0.1849955 1 0.0421696503 0.6011656 0.1949709 1 0.1000000000 0.6020123 0.1966378 3 0.0000000000 0.5994884 0.1909870 3 0.0001000000 0.5994841 0.1898289 3 0.0002371374 0.6026696 0.1972467 3 0.0005623413 0.6023385 0.1968234 3 0.0013335214 0.6014989 0.1952293 3 0.0031622777 0.5930657 0.1797088 3 0.0074989421 0.6094037 0.2115621 3 0.0177827941 0.6053659 0.2038492 3 0.0421696503 0.6063616 0.2069790 3 0.1000000000 0.6062154 0.2053208 5 0.0000000000 0.5954286 0.1824067 5 0.0001000000 0.5999954 0.1920286 5 0.0002371374 0.5929118 0.1784488 5 0.0005623413 0.5942436 0.1811213 5 0.0013335214 0.5976199 0.1879154 5 0.0031622777 0.5998076 0.1921799 5 0.0074989421 0.6023336 0.1972625 5 0.0177827941 0.6058581 0.2043481 5 0.0421696503 0.6107497 0.2151394 5 0.1000000000 0.6053736 0.2041157 7 0.0000000000 0.5792771 0.1534069 7 0.0001000000 0.5868479 0.1674675 7 0.0002371374 0.5922169 0.1779055 7 0.0005623413 0.5866473 0.1664921 7 0.0013335214 0.5971093 0.1870028 7 0.0031622777 0.5908964 0.1745331 7 0.0074989421 0.6055304 0.2042054 7 0.0177827941 0.6092053 0.2115698 7 0.0421696503 0.6082266 0.2104728 7 0.1000000000 0.6040298 0.2010179 9 0.0000000000 0.5871833 0.1693874 9 0.0001000000 0.5875258 0.1681067 9 0.0002371374 0.5844843 0.1632334 9 0.0005623413 0.5822853 0.1579851 9 0.0013335214 0.5827973 0.1586526 9 0.0031622777 0.5873324 0.1679448 9 0.0074989421 0.5901888 0.1732924 9 0.0177827941 0.6043405 0.2013014 9 0.0421696503 0.6122553 0.2182269 9 0.1000000000 0.6035178 0.1999213 11 0.0000000000 0.5794367 0.1522612 11 0.0001000000 0.5853317 0.1644861 11 0.0002371374 0.5836454 0.1624478 11 0.0005623413 0.5865146 0.1678897 11 0.0013335214 0.5861693 0.1656024 11 0.0031622777 0.5908936 0.1747307 11 0.0074989421 0.5950925 0.1833357 11 0.0177827941 0.6043313 0.2016203 11 0.0421696503 0.6105817 0.2149566 11 0.1000000000 0.6031817 0.1992942 13 0.0000000000 0.5888767 0.1714118 13 0.0001000000 0.5843403 0.1620815 13 0.0002371374 0.5864998 0.1665880 13 0.0005623413 0.5829696 0.1591928 13 0.0013335214 0.5821321 0.1582996 13 0.0031622777 0.5930608 0.1796795 13 0.0074989421 0.5935700 0.1802362 13 0.0177827941 0.6021641 0.1968581 13 0.0421696503 0.6107476 0.2154965 13 0.1000000000 0.6046957 0.2025467 15 0.0000000000 0.5818080 0.1560407 15 0.0001000000 0.5770560 0.1491584 15 0.0002371374 0.5763951 0.1467020 15 0.0005623413 0.5826101 0.1593395 15 0.0013335214 0.5900391 0.1734328 15 0.0031622777 0.5851608 0.1633889 15 0.0074989421 0.5944159 0.1818350 15 0.0177827941 0.6039994 0.2004030 15 0.0421696503 0.6100704 0.2141646 15 0.1000000000 0.6043560 0.2020346 17 0.0000000000 0.5801111 0.1553176 17 0.0001000000 0.5748789 0.1434048 17 0.0002371374 0.5740541 0.1426087 17 0.0005623413 0.5815868 0.1578588 17 0.0013335214 0.5794387 0.1531597 17 0.0031622777 0.5823235 0.1579740 17 0.0074989421 0.5878351 0.1688083 17 0.0177827941 0.6041753 0.2011916 17 0.0421696503 0.6100767 0.2140261 17 0.1000000000 0.6058707 0.2049583 19 0.0000000000 0.5769198 0.1483938 19 0.0001000000 0.5818002 0.1581756 19 0.0002371374 0.5797791 0.1543552 19 0.0005623413 0.5785913 0.1513950 19 0.0013335214 0.5804387 0.1545173 19 0.0031622777 0.5854843 0.1644059 19 0.0074989421 0.5934054 0.1799900 19 0.0177827941 0.6023216 0.1973172 19 0.0421696503 0.6105739 0.2151288 19 0.1000000000 0.6033462 0.2000386 Accuracy was used to select the optimal model using the largest value. The final values used for the model were size = 9 and decay = 0.04216965.
Visualize ANN [this doesn't work :(]
library(caret) source_url('https://gist.githubusercontent.com/fawda123/7471137/raw/466c1474d0a505ff044412703516c34f1a4684a5/nnet_plot_update.r') plot.nnet(fn.nn)
Predict test data
fn.nn.predict <- predict(fn.nn, newdata = x.test) confusionMatrix(fn.nn.predict,as.factor(y.test))
Confusion Matrix and Statistics Reference Prediction DOWN UP DOWN 95 72 UP 23 51 Accuracy : 0.6058 95% CI : (0.541, 0.6679) No Information Rate : 0.5104 P-Value NIR] : 0.001806 Kappa : 0.2179 Mcnemar's Test P-Value : 8.449e-07 Sensitivity : 0.8051 Specificity : 0.4146 Pos Pred Value : 0.5689 Neg Pred Value : 0.6892 Prevalence : 0.4896 Detection Rate : 0.3942 Detection Prevalence : 0.6929 Balanced Accuracy : 0.6099 'Positive' Class : DOWN
8.3.2 ANN for Regression
Pre-process data to use fnx classification problem dataset split into test and train
xy <- fnx xy$Expected <- lead(xy$FTSE.Close) xy<- xy[-nrow(xy),] xy.train <- head(xy, -60) xy.test <- tail(xy, 60)
Train on training set
set.seed(12345) ctrl <- trainControl(method = "repeatedcv", repeats = 3, number = 5) fnReg.nn <- train(Expected ~ ., data = xy.train, method = "nnet", trControl = ctrl, tuneLength = 5, trace = FALSE) fnReg.nn
Neural Network 1368 samples 10 predictor No pre-processing Resampling: Cross-Validated (5 fold, repeated 3 times) Summary of sample sizes: 1094, 1094, 1094, 1095, 1095, 1093, ... Resampling results across tuning parameters: size decay RMSE Rsquared MAE 1 0e+00 0.03001576 0.9847348 0.02276178 1 1e-04 0.03013584 0.9846119 0.02285315 1 1e-03 0.03095296 0.9838156 0.02354345 1 1e-02 0.03367255 0.9811078 0.02573516 1 1e-01 0.04162665 0.9742108 0.03090446 3 0e+00 0.02550131 0.9888659 0.01880705 3 1e-04 0.02576956 0.9886560 0.01926426 3 1e-03 0.02543368 0.9889975 0.01895336 3 1e-02 0.02957528 0.9853183 0.02244243 3 1e-01 0.03746831 0.9782222 0.02789070 5 0e+00 0.02506330 0.9892662 0.01846131 5 1e-04 0.02523152 0.9891421 0.01871194 5 1e-03 0.02558013 0.9888594 0.01894277 5 1e-02 0.02928744 0.9855734 0.02225197 5 1e-01 0.03625234 0.9793464 0.02694509 7 0e+00 0.02513850 0.9892017 0.01852413 7 1e-04 0.02512707 0.9892271 0.01853013 7 1e-03 0.02551769 0.9889172 0.01898627 7 1e-02 0.02924782 0.9856128 0.02224497 7 1e-01 0.03586885 0.9796820 0.02663118 9 0e+00 0.02511433 0.9892187 0.01848913 9 1e-04 0.02511768 0.9892296 0.01855158 9 1e-03 0.02551788 0.9889041 0.01895170 9 1e-02 0.02916250 0.9856959 0.02218335 9 1e-01 0.03543371 0.9801297 0.02631143 RMSE was used to select the optimal model using the smallest value. The final values used for the model were size = 5 and decay = 0.
Visualize [doesn't work]
plot.nnet(fnReg.nn)
Predict test set
fnReg.predict <- predict(fnReg.nn, newdata = xy.test) modelValues <- data.frame(obs = xy.test$Expected, pred = fnReg.predict) defaultSummary(modelValues)
Error in predict(fnReg.nn, newdata = xy.test) : object 'fnReg.nn' not found RMSE Rsquared MAE 6262.2829355 0.8099375 6260.6262660
Plot (normalized)
plot(xy.test$Expected, type = "l", ylim = c(0, 1), ylab = "Normalized FTSE Price") lines(fnReg.predict, col = "red") legend("topright", c("Test Data", "Prediction"), col = c("black", "red"), lty = c(1, 1))
Plot (original values)
FTSE.min <- min(fn$FTSE.Close) FTSE.max <- max(fn$FTSE.Close) expected <- FTSE.min + xy.test$Expected*(FTSE.max - FTSE.min) predicted <- FTSE.min + fnReg.predict*(FTSE.max - FTSE.min) plot(expected, type = "l", ylim = c(4500, 7500), ylab = "FTSE Price") lines(predicted, col = "red") #legend("topright", c("Expected", "Predicted"), col = c("black", "red"), #lty = c(1, 1))
Calculate RMSE and R-squared for predicted price for original data
modelValues <- data.frame(obs = expected, pred = predicted) defaultSummary(modelValues)
RMSE Rsquared MAE 1.439217e+07 8.099375e-01 1.438836e+07
9 Trading Strategies and Backtesting
9.1 Trading Strategy Identification
A mean-reverting strategy exploits the fact that price returns to its long-term average price after a temporary deviation: buy low and sell high Momentum starategies based on strong moves int he market, where more direction follows direction: higher highs and lower lows Trend-trading stratggies buy stocks on upward trend and sell on downward trend.
To illustrate:
library(quantmod) getSymbols("IBM", from = "2012-01-01", to = "2013-08-01") chartSeries(IBM, theme = "white") #, minor.ticks = F) addTA(SMA(Cl(IBM), n = 40), on = 1, col = "blue")
Another example using GS
getSymbols("GS", from = "2011-01-01", to = "2015-07-31") gs <- GS["2012-07-31::2015-07-31"] chartSeries(gs, theme = "white")#, minor.ticks = F) addTA(SMA(Cl(GS), n = 200), on = 1, col = "blue")
9.2 Trading Signals
9.2.1 Signals from moving average
When the price deviates more than SignalIn standard deviations from mean, we enter position. SignalIn is a free parameter that can be optimized during training. The length of the moving window (for mean and standard deviation) can also be optimized. We exit positions when price returns to SignalOut standard deviations from the mean. SignalOut is usually set to 0 (the mean). When moving window is short, SignalIn and SignalOut will be small, and we will have shorter holding periods and thus more trades.
Bollinger bands IBM
chartSeries(IBM, theme = "white", yrange = c(170, 230), TA = "addBBands(n = 40, maType = 'SMA')")
Strategy Helper file
Calculate moving average z-scores using strategyHelper maZscore function
zscore <- maZscore(Cl(IBM), 40) plot(zscore)
9.2.2 Signals from Linear Regression
This is price - pred.price / stdev This is also a rolling signal.
[error: invalid 'y' type in 'x & & y']
source("./R/strategyHelper.R") zscore <- lrZscore(as.numeric(Cl(IBM)), as.numeric(40)) plot(zscore$Zscore, main = "Linear regression Zscore for IBM") plot(Cl(IBM)) lines(zscore$FittedValues, col="red")
9.2.3 Signals from RSI
Relative Strength Index (RSI) Is a momentum oscillator that measures spped and change of price movement, 0 to 100. Considered overbought when > 70 and oversold when below 30. The RSI formulas = 100 - 100 / (1 + RS'), where RS = Average Gain / Average Loss Losses are expressed as positive values. To extract trading signal, re-normalize so it's bettwenn -2 and 2 Signal = RSI - 50 / 25
Apply RSI signal
sig <- rsiSignal(Cl(IBM), 14) plot(sig)
9.2.4 Signal from Williams %R
Williams %R is a momentum indicator that show the level of hte close price relative to the ihghest high for a given moving window. It oscillates from -100 to 0. Readings from 0 to -20 are overbough, while reading from -80 to -100 are considered oversold. The formula is WilliamsR = -100 * (highest high - Close) / (highest high - lowest low) Trading signal after renormalize to be between -2 and 2: Signal = WilliamsR/25+2 = -4 * (highest high - Close) / (highest high - lowest low) + 2
Applying Williams %R (moving window of 14)
wpr <- wprSignal(IBM[, c("IBM.High", "IBM.Low", "IBM.Close")], 14) plot(wpr)
9.3 Backtest System Implementation
Use trading signal to test strategies on historical data. Can calculate daily and accumulated profit and loss and other risk measures for a given trading strategy. Assumptions here: no stop-loss, ignore transaction costs, ignore bid/ask spreads, and assumme ordered can be filled fully at specified close price. We have to use the popular TWAP and or VWAP to execute our orders.
9.3.1 Profit and Loss Computation
The first long-short trading strategy is located in strategyHelper.R This strategies takes as its input price, zin, and zout. It takes a long position (+1) when zscore is less than -zin and flat the position when zscore is greater than -zout. It takes a short position (-1) when zscore is greater than zin and flat the position when zscore < zout. It exits either position when zout = 0. The first day is flat (0) and the last day is flat (0). This function also replaces each na with the previous non-na value.
There are columns for Long position and columns for Shrot position.
The position is calculated as the lag of the sum of Long and Short positions. This lag is important to avoid look-forward bias.
This also calculates Equity, REturn, ReturnDaily, Return Hold, and Return Daily Hold.
Test the longshort Trading strategy based on moving average
source("./R/strategyHelper.R") zscore <- maZscore(Cl(IBM), 40) price <- na.omit(data.frame(Price = Cl(IBM), Zscore <- zscore[,1])) strategy <- longShortStrategy(price, 2, 0) tail(strategy)
Error in diff(strategy$Position, na.pad = FALSE) : object 'strategy' not found Error in tail(strategy) : object 'strategy' not found
Plot cumulative returns
plot(strategy$Return, main = "Cumulative returns for IBM", ylim = c(-0.1, 0.27)) lines(strategy$RetHold, col = "red")
source("./R/strategyHelper.R") zscore <- lrZscore(Cl(IBM), 40) price <- na.omit(data.frame(Price = Cl(IBM), Zscore <- zscore[,1])) strategy <- longShortStrategy(price, 2, 0) tail(strategy)
Error in by.column && !is.null(dim(data)) : invalid 'x' type in 'x && y' Error in diff(strategy$Position, na.pad = FALSE) : object 'strategy' not found Error in tail(strategy) : object 'strategy' not found
plot(strategy$Return, main = "Cumulative returns for IBM", ylim = c(-0.1, 0.27)) lines(strategy$RetHold, col = "red")
calculate other metrics if interested
numTrades <- sum(abs(diff(strategy$Position, na.pad = FALSE))) numtrades Pnl.oneShare <- sum(diff(strategy$Price) * strategy$Position, na.rm = TRUE) Pnl.oneShare #P&L for one share of investment rtn <- 100 * sum(ROC(strategy$Price) * strategy$Position, na.rm=TRUE) rtn #Percentage return of fully invested account
Error in diff(strategy$Position, na.pad = FALSE) : object 'strategy' not found Error: object 'numtrades' not found Error in diff(strategy$Price) : object 'strategy' not found Error: object 'Pnl.oneShare' not found Error in is.xts(x) : object 'strategy' not found GSPC.Close 2012-06-01 NA 2012-06-04 1.095485e-04 2012-06-05 5.710514e-03 2012-06-06 2.278778e-02 2012-06-07 -1.064704e-04 2012-06-08 8.081421e-03 2012-06-11 -1.270042e-02 2012-06-12 1.158339e-02 2012-06-13 -7.048030e-03 2012-06-14 1.075659e-02 2012-06-15 1.028475e-02 2012-06-18 1.443704e-03 2012-06-19 9.767834e-03 2012-06-20 -1.687780e-03 2012-06-21 -2.251321e-02 2012-06-22 7.148989e-03 2012-06-25 -1.608350e-02 2012-06-26 4.761368e-03 2012-06-27 8.944782e-03 2012-06-28 -2.112029e-03 2012-06-29 2.461479e-02 2012-07-02 2.456292e-03 2012-07-03 6.212772e-03 2012-07-05 -4.698042e-03 2012-07-06 -9.477418e-03 2012-07-09 -1.640176e-03 2012-07-10 -8.159122e-03 2012-07-11 -1.492404e-05 2012-07-12 -4.999574e-03 2012-07-13 1.636276e-02 2012-07-16 -2.316995e-03 2012-07-17 7.382356e-03 2012-07-18 6.658275e-03 2012-07-19 2.713416e-03 2012-07-20 -1.011262e-02 2012-07-23 -8.948979e-03 2012-07-24 -9.082050e-03 2012-07-25 -3.139108e-04 2012-07-26 1.640566e-02 2012-07-27 1.890081e-02 2012-07-30 -4.834765e-04 2012-07-31 -4.326173e-03 2012-08-01 -2.904193e-03 2012-08-02 -7.531963e-03 2012-08-03 1.886129e-02 2012-08-06 2.326561e-03 2012-08-07 5.093763e-03 2012-08-08 6.206337e-04 2012-08-09 4.135999e-04 2012-08-10 2.186050e-03 2012-08-13 -1.252685e-03 2012-08-14 -1.281542e-04 2012-08-15 1.138991e-03 2012-08-16 7.075421e-03 2012-08-17 1.870383e-03 2012-08-20 -2.117484e-05 2012-08-21 -3.503667e-03 2012-08-22 2.263774e-04 2012-08-23 -8.104999e-03 2012-08-24 6.433988e-03 2012-08-27 -4.891348e-04 2012-08-28 -8.085086e-04 2012-08-29 8.439927e-04 2012-08-30 -7.836430e-03 2012-08-31 5.060470e-03 2012-09-04 -1.166640e-03 2012-09-05 -1.068232e-03 2012-09-06 2.022954e-02 2012-09-07 4.041795e-03 2012-09-10 -6.166806e-03 2012-09-11 3.130052e-03 2012-09-12 2.090506e-03 2012-09-13 1.617817e-02 2012-09-14 3.951136e-03 2012-09-17 -3.129583e-03 2012-09-18 -1.280595e-03 2012-09-19 1.184852e-03 2012-09-20 -5.408799e-04 2012-09-21 -7.532230e-05 2012-09-24 -2.235150e-03 2012-09-25 -1.055739e-02 2012-09-26 -5.753254e-03 2012-09-27 9.602728e-03 2012-09-28 -4.487808e-03 2012-10-01 2.647997e-03 2012-10-02 8.719069e-04 2012-10-03 3.617857e-03 2012-10-04 7.148822e-03 2012-10-05 -3.216406e-04 2012-10-08 -3.462724e-03 2012-10-09 -9.940183e-03 2012-10-10 -6.207255e-03 2012-10-11 1.953703e-04 2012-10-12 -2.970545e-03 2012-10-15 8.045471e-03 2012-10-16 1.021756e-02 2012-10-17 4.108606e-03 2012-10-18 -2.446720e-03 2012-10-19 -1.671015e-02 2012-10-22 4.394858e-04 2012-10-23 -1.454924e-02 2012-10-24 -3.090152e-03 2012-10-25 2.991065e-03 2012-10-26 -7.292481e-04 2012-10-31 1.558677e-04 2012-11-01 1.086721e-02 2012-11-02 -9.423718e-03 2012-11-05 2.161472e-03 2012-11-06 7.822509e-03 2012-11-07 -2.399048e-02 2012-11-08 -1.227993e-02 2012-11-09 1.697251e-03 2012-11-12 1.304789e-04 2012-11-13 -3.993383e-03 2012-11-14 -1.394887e-02 2012-11-15 -1.594816e-03 2012-11-16 4.828275e-03 2012-11-19 1.966738e-02 2012-11-20 6.631665e-04 2012-11-21 2.317493e-03 2012-11-23 1.294220e-02 2012-11-26 -2.031644e-03 2012-11-27 -5.240294e-03 2012-11-28 7.825331e-03 2012-11-29 4.260553e-03 2012-11-30 1.624947e-04 2012-12-03 -4.756519e-03 2012-12-04 -1.711276e-03 2012-12-05 1.583607e-03 2012-12-06 3.301136e-03 2012-12-07 2.916662e-03 2012-12-10 3.385036e-04 2012-12-11 6.527531e-03 2012-12-12 4.481389e-04 2012-12-13 -6.341490e-03 2012-12-14 -4.143976e-03 2012-12-17 1.180069e-02 2012-12-18 1.142119e-02 2012-12-19 -7.618145e-03 2012-12-20 5.473104e-03 2012-12-21 -9.422944e-03 2012-12-24 -2.443279e-03 2012-12-26 -4.798957e-03 2012-12-27 -1.219185e-03 2012-12-28 -1.111145e-02 2012-12-31 1.680003e-02 2013-01-02 2.508612e-02 2013-01-03 -2.087796e-03 2013-01-04 4.853300e-03 2013-01-07 -3.128003e-03 2013-01-08 -3.247640e-03 2013-01-09 2.652346e-03 2013-01-10 7.568700e-03 2013-01-11 -4.751492e-05 2013-01-14 -9.311048e-04 2013-01-15 1.128033e-03 2013-01-16 1.969725e-04 2013-01-17 5.627060e-03 2013-01-18 3.397492e-03 2013-01-22 4.418332e-03 2013-01-23 1.506342e-03 2013-01-24 6.614196e-06 2013-01-25 5.430709e-03 2013-01-28 -1.851334e-03 2013-01-29 5.093004e-03 2013-01-30 -3.907245e-03 2013-01-31 -2.566592e-03 2013-02-01 1.000251e-02 2013-02-04 -1.160583e-02 2013-02-05 1.036263e-02 2013-02-06 5.490198e-04 2013-02-07 -1.807031e-03 2013-02-08 5.641995e-03 2013-02-11 -6.063013e-04 2013-02-12 1.594001e-03 2013-02-13 5.920875e-04 2013-02-14 6.904334e-04 2013-02-15 -1.045628e-03 2013-02-19 7.309694e-03 2013-02-20 -1.248171e-02 2013-02-21 -6.323005e-03 2013-02-22 8.734214e-03 2013-02-25 -1.847928e-02 2013-02-26 6.090876e-03 2013-02-27 1.264570e-02 2013-02-28 -8.644531e-04 2013-03-01 2.321159e-03 2013-03-04 4.600127e-03 2013-03-05 9.520552e-03 2013-03-06 1.083925e-03 2013-03-07 1.814844e-03 2013-03-08 4.471129e-03 2013-03-11 3.243819e-03 2013-03-12 -2.406146e-03 2013-03-13 1.313190e-03 2013-03-14 5.587352e-03 2013-03-15 -1.619774e-03 2013-03-18 -5.525570e-03 2013-03-19 -2.425470e-03 2013-03-20 6.675164e-03 2013-03-21 -8.316924e-03 2013-03-22 7.148644e-03 2013-03-25 -3.345630e-03 2013-03-26 7.754964e-03 2013-03-27 -5.885231e-04 2013-03-28 4.048463e-03 2013-04-01 -4.483617e-03 2013-04-02 5.158934e-03 2013-04-03 -1.060213e-02 2013-04-04 4.040279e-03 2013-04-05 -4.304145e-03 2013-04-08 6.282959e-03 2013-04-09 3.538066e-03 2013-04-10 1.211544e-02 2013-04-11 3.545956e-03 2013-04-12 -2.840798e-03 2013-04-15 -2.323413e-02 2013-04-16 1.420584e-02 2013-04-17 -1.443131e-02 2013-04-18 -6.723556e-03 2013-04-19 8.808989e-03 2013-04-22 4.650798e-03 2013-04-23 1.036531e-02 2013-04-24 6.340319e-06 2013-04-25 4.026614e-03 2013-04-26 -1.843812e-03 2013-04-29 7.160315e-03 2013-04-30 2.481817e-03 2013-05-01 -9.351473e-03 2013-05-02 9.364004e-03 2013-05-03 1.047956e-02 2013-05-06 1.905961e-03 2013-05-07 5.216639e-03 2013-05-08 4.130539e-03 2013-05-09 -3.693918e-03 2013-05-10 4.312344e-03 2013-05-13 4.288884e-05 2013-05-14 1.009107e-02 2013-05-15 5.101103e-03 2013-05-16 -5.022332e-03 2013-05-17 1.024741e-02 2013-05-20 -7.078686e-04 2013-05-21 1.720905e-03 2013-05-22 -8.308074e-03 2013-05-23 -2.928115e-03 2013-05-24 -5.515174e-04 2013-05-28 6.320962e-03 2013-05-29 -7.072937e-03 2013-05-30 3.663625e-03 2013-05-31 -1.441058e-02 2013-06-03 5.918440e-03 2013-06-04 -5.526048e-03 2013-06-05 -1.387555e-02 2013-06-06 8.454455e-03 2013-06-07 1.274991e-02 2013-06-10 -3.468734e-04 2013-06-11 -1.020526e-02 2013-06-12 -8.404777e-03 2013-06-13 1.467607e-02 2013-06-14 -5.902401e-03 2013-06-17 7.538876e-03 2013-06-18 7.760964e-03 2013-06-19 -1.394830e-02 2013-06-20 -2.532842e-02 2013-06-21 2.666219e-03 2013-06-24 -1.221937e-02 2013-06-25 9.452456e-03 2013-06-26 9.544790e-03 2013-06-27 6.180691e-03 2013-06-28 -4.298789e-03 2013-07-01 5.389200e-03 2013-07-02 -5.450568e-04 2013-07-03 8.237078e-04 2013-07-05 1.015005e-02 2013-07-08 5.237805e-03 2013-07-09 7.203662e-03 2013-07-10 1.815761e-04 2013-07-11 1.346321e-02 2013-07-12 3.081730e-03 2013-07-15 1.373935e-03 2013-07-16 -3.715655e-03 2013-07-17 2.770206e-03 2013-07-18 5.020342e-03 2013-07-19 1.608756e-03 2013-07-22 2.030962e-03 2013-07-23 -1.853653e-03 2013-07-24 -3.818503e-03 2013-07-25 2.553210e-03 2013-07-26 8.279512e-04 2013-07-29 -3.743034e-03 2013-07-30 3.737472e-04 2013-07-31 -1.364188e-04 2013-08-01 1.246259e-02 2013-08-02 1.639114e-03 2013-08-05 -1.480931e-03 2013-08-06 -5.739473e-03 2013-08-07 -3.813126e-03 2013-08-08 3.877921e-03 2013-08-09 -3.576348e-03 2013-08-12 -1.153586e-03 2013-08-13 2.772210e-03 2013-08-14 -5.190063e-03 2013-08-15 -1.438457e-02 2013-08-16 -3.310067e-03 2013-08-19 -5.917778e-03 2013-08-20 3.813913e-03 2013-08-21 -5.796370e-03 2013-08-22 8.582442e-03 2013-08-23 3.939242e-03 2013-08-26 -4.047839e-03 2013-08-27 -1.600154e-02 2013-08-28 2.743878e-03 2013-08-29 1.961477e-03 2013-08-30 -3.179367e-03 2013-09-03 4.155575e-03 2013-09-04 8.084188e-03 2013-09-05 1.209132e-03 2013-09-06 5.442973e-05 2013-09-09 9.943282e-03 2013-09-10 7.318940e-03 2013-09-11 3.047635e-03 2013-09-12 -3.386142e-03 2013-09-13 2.711001e-03 2013-09-16 5.677009e-03 2013-09-17 4.208869e-03 2013-09-18 1.210412e-02 2013-09-19 -1.844654e-03 2013-09-20 -7.243054e-03 2013-09-23 -4.730759e-03 2013-09-24 -2.600522e-03 2013-09-25 -2.743225e-03 2013-09-26 3.479366e-03 2013-09-27 -4.082122e-03 2013-09-30 -6.047480e-03 2013-10-01 7.966725e-03 2013-10-02 -6.668919e-04 2013-10-03 -9.019973e-03 2013-10-04 7.028467e-03 2013-10-07 -8.542748e-03 2013-10-08 -1.240875e-02 2013-10-09 5.737416e-04 2013-10-10 2.159562e-02 2013-10-11 6.266595e-03 2013-10-14 4.066441e-03 2013-10-15 -7.088790e-03 2013-10-16 1.373280e-02 2013-10-17 6.721314e-03 2013-10-18 6.527406e-03 2013-10-21 9.173211e-05 2013-10-22 5.721116e-03 2013-10-23 -4.735753e-03 2013-10-24 3.252838e-03 2013-10-25 4.385215e-03 2013-10-28 1.328816e-03 2013-10-29 5.568661e-03 2013-10-30 -4.887849e-03 2013-10-31 -3.846771e-03 2013-11-01 2.899215e-03 2013-11-04 3.564199e-03 2013-11-05 -2.809531e-03 2013-11-06 4.256468e-03 2013-11-07 -1.327044e-02 2013-11-08 1.333821e-02 2013-11-11 7.226705e-04 2013-11-12 -2.373206e-03 2013-11-13 8.062753e-03 2013-11-14 4.825597e-03 2013-11-15 4.213147e-03 2013-11-18 -3.705053e-03 2013-11-19 -2.045055e-03 2013-11-20 -3.642236e-03 2013-11-21 8.095706e-03 2013-11-22 4.949190e-03 2013-11-25 -1.264141e-03 2013-11-26 1.497935e-04 2013-11-27 2.481998e-03 2013-11-29 -7.859980e-04 2013-12-02 -2.722724e-03 2013-12-03 -3.197956e-03 2013-12-04 -1.304343e-03 2013-12-05 -4.349016e-03 2013-12-06 1.117520e-02 2013-12-09 1.815451e-03 2013-12-10 -3.184725e-03 2013-12-11 -1.138140e-02 2013-12-12 -3.777689e-03 2013-12-13 -1.014154e-04 2013-12-16 6.300152e-03 2013-12-17 -3.105806e-03 2013-12-18 1.651091e-02 2013-12-19 -5.800970e-04 2013-12-20 4.807155e-03 2013-12-23 5.304029e-03 2013-12-24 2.911504e-03 2013-12-26 4.734305e-03 2013-12-27 -3.366415e-04 2013-12-30 -1.792699e-04 2013-12-31 3.951856e-03 2014-01-02 -8.901413e-03 2014-01-03 -3.330203e-04 2014-01-06 -2.514927e-03 2014-01-07 6.063345e-03 2014-01-08 -2.122317e-04 2014-01-09 3.482487e-04 2014-01-10 2.304030e-03 2014-01-13 -1.265597e-02 2014-01-14 1.075988e-02 2014-01-15 5.152889e-03 2014-01-16 -1.348028e-03 2014-01-17 -3.902781e-03 2014-01-21 2.769912e-03 2014-01-22 5.746998e-04 2014-01-23 -8.929325e-03 2014-01-24 -2.109642e-02 2014-01-27 -4.888222e-03 2014-01-28 6.121875e-03 2014-01-29 -1.026170e-02 2014-01-30 1.120404e-02 2014-01-31 -6.486290e-03 2014-02-03 -2.309660e-02 2014-02-04 7.612043e-03 2014-02-05 -2.030282e-03 2014-02-06 1.236305e-02 2014-02-07 1.321419e-02 2014-02-10 1.568005e-03 2014-02-11 1.100138e-02 2014-02-12 -2.692985e-04 2014-02-13 5.793212e-03 2014-02-14 4.797689e-03 2014-02-18 1.157803e-03 2014-02-19 -6.545862e-03 2014-02-20 6.013342e-03 2014-02-21 -1.920567e-03 2014-02-24 6.167455e-03 2014-02-25 -1.348591e-03 2014-02-26 2.169971e-05 2014-02-27 4.935881e-03 2014-02-28 2.778824e-03 2014-03-03 -7.405866e-03 2014-03-04 1.515232e-02 2014-03-05 -5.335244e-05 2014-03-06 1.716933e-03 2014-03-07 5.379447e-04 2014-03-10 -4.633536e-04 2014-03-11 -5.095097e-03 2014-03-12 3.051242e-04 2014-03-13 -1.177009e-02 2014-03-14 -2.825766e-03 2014-03-17 9.567718e-03 2014-03-18 7.193682e-03 2014-03-19 -6.150525e-03 2014-03-20 6.022334e-03 2014-03-21 -2.936980e-03 2014-03-24 -4.876581e-03 2014-03-25 4.394271e-03 2014-03-26 -7.024937e-03 2014-03-27 -1.901892e-03 2014-03-28 4.629489e-03 2014-03-31 7.892872e-03 2014-04-01 7.014688e-03 2014-04-02 2.849263e-03 2014-04-03 -1.127085e-03 2014-04-04 -1.261654e-02 2014-04-07 -1.080831e-02 2014-04-08 3.743538e-03 2014-04-09 1.085904e-02 2014-04-10 -2.110597e-02 2014-04-11 -9.532060e-03 2014-04-14 8.183708e-03 2014-04-15 6.734579e-03 2014-04-16 1.043387e-02 2014-04-17 1.362924e-03 2014-04-21 3.768016e-03 2014-04-22 4.083789e-03 2014-04-23 -2.215767e-03 2014-04-24 1.715488e-03 2014-04-25 -8.129345e-03 2014-04-28 3.230812e-03 2014-04-29 4.749460e-03 2014-04-30 2.987550e-03 2014-05-01 -1.432715e-04 2014-05-02 -1.349355e-03 2014-05-05 1.869467e-03 2014-05-06 -9.029031e-03 2014-05-07 5.600755e-03 2014-05-08 -1.374569e-03 2014-05-09 1.518323e-03 2014-05-12 9.626256e-03 2014-05-13 4.216689e-04 2014-05-14 -4.712090e-03 2014-05-15 -9.405904e-03 2014-05-16 3.739962e-03 2014-05-19 3.837414e-03 2014-05-20 -6.519605e-03 2014-05-21 8.083340e-03 2014-05-22 2.359444e-03 2014-05-23 4.239393e-03 2014-05-27 5.969950e-03 2014-05-28 -1.114693e-03 2014-05-29 5.352759e-03 2014-05-30 1.841980e-03